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Models  of  The  Near-Space 
Geophysical  Environment 


1.  INTRODUCTION 


This  report  provides  descriptions  of  environmental  models  used  in  the  near-earth  environment  for  the 
analysis  of  earth-based  and  satellite  data.  Included  are  models  of  the  ionosphere,  atmosphere, 
radiation  belt.  Earth’s  magnetic  field,  and  orbital  prediction.  These  models  are  used  at  PL  for  a 
variety  of  purposes,  in  particular,  for  the  analysis  of  experimental  data  obtained  fi'om  ground  based 
instruments  and,  where  applicable,  instruments  carried  on  balloons,  rockets,  and  satellites.  Some  of 
the  models  described  here  were  developed  at  PL  based  upon  such  measurements.  Extensive 
literature  exists  describing  the  scientific  background  of  the  models.  The  Handbook  of  Geophysics 
and  the  Space  Environment  [Jursa,  1985]  is  in  particular  pertinent  to  PL  and  this  report,  but  clearly 
shows  that  many  models  remain  to  be  documented.  For  many  of  the  models  described  here,  computer 
output  is  provided  together  with  estimates  of  computer  (CPU)  times  required  to  generate  such 
output.  The  computer  codes  (subroutines)  for  these  models  are  frequently  incorporated  into  programs 
used  for  data  analysis  and  reduction,  as  well  as  mission  planning.  It  is  hoped  that  this  report  will  serve 
as  a  useful  reference  for  planning  such  future  missions.  Although  most  of  the  programs  are  in  the 
public  domain,  some  of  the  ionospheric  and  radiation  belt  models  are  developed  and  distributed  by 
the  PL  research  branches. 
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2.  IONOSPHERE  MODELS  AND  IONOSPHERE  RAY  TRACING 


This  section  provides  a  description  of  several  ionospheric  models,  which  provide  empirical  and 
mathematical  representations  of  the  electron  density  of  the  ionosphere,  and,  in  some  cases,  density 
of  ions,  electron  and  ion  temperatures  for  a  variety  of  conditions  ranging  from  magnetically  quiet, 
to  active  conditions.  Some  of  the  models  provide  estimates  of  ionospheric  propagation  conditions 
based  upon  solar  activity,  which  can  be  used  as  the  basis  of  ionospheric  ray  tracing  routines. 
Descriptions  of  two  ray  tracing  routines  are  also  provided. 


2.1  International  Reference  Ionosphere  1990  (IRI-90) 

The  International  Reference  Ionosphere  (IRI)  [Bilitza,  1990]  is  the  empirical  reference  model  of 
ionospheric  densities  and  temperatures  (electrons  and  ions)  recommended  for  international  use  by  the 
Committee  on  Space  Research  (COSPAR)  and  the  International  Union  of  Radio  Science  (URSI). 
IRI  was  established  in  a  world-wide  collaborative  effort  beginning  in  the  late  sixties,  and  is  updated 
bi-yearly  during  special  workshops. 

IRI  is  based  on  both  ground-based  (ionosonde,  absorption,  and  incoherent  scatter)  and  spacecraft 
measurements  (Alouette,  ISIS,  AE,  AEROS,  DE,  and  rockets).  The  model  provides  monthly  mean 
values  for  magnetically  quiet  conditions  at  non-auroral  latitudes,  for  altitudes  ranging  from  50  km  to 
2000  km.  This  model  serves  as  a  standard  reference  for  various  purposes,  such  as  design  of 
experiments,  estimation  of  environmental  and  other  effects,  and  testing  theories.  IRI  calculates 
electron  density,  electron  and  ion  temperatures,  and  the  composition  of  positive  ions.  It  computes 
the  density  for  atomic  ions  O'",  H'",  He"",  and  for  molecular  O2  and  NO'". 

IRI-90  offers  several  new  options  not  found  in  its  predecessor,  IRI-86: 

•  URSI-1989  foF2  model 

•  Gulyaeva- 1987  model  for  F2  bottomside  thickness 

•  Analytical  LAY-representation  of  E-F  region 

•  Danilov-Yaichnikov-1985  model  for  ion  composition 

A  significant  improvement  in  the  representation  of  the  foE  nighttime  variation  has  been  incorporated 
in  IRI-90.  In  addition,  the  neutral  temperature  is  obtained  with  CIRA-86,  rather  than  with  CIRA-72. 

IRI  will  produce  profiles  in  latitude  or  longitude  (geodetic  or  geomagnetic),  solar  activity,  month,  day 
of  year,  local  (or  universal)  time,  or  altitude.  The  user  must  supply  latitude,  longitude,  sunspot 
number,  and  time.  Calculation  of  the  full  set  of  output  parameters  for  a  profile  containing  9000  points 
requires  20  seconds  on  a  66  MHz  486  PC.  Table  1  shows  a  sample  IRI  profile  for  12  UT  on  October 
15,  at  42.5°N  and  -71.3®E,  with  a  sunspot  number  of  100. 


Table  1 .  Sample  IRI-90  Output 


H(KM) 

ELECTRON  DENSITY 

TEMPERATURES 

ION 

PERCENTAGE 

DENSITIES 

NE/CM-3 

NE/NMF2 

TN(K) 

TI  (K) 

TE(K) 

TE/TI 

0+ 

H+ 

He+ 

02+ 

NO+ 

100.0 

53912 

.0936 

-1 

-1 

-1  • 

-1.00 

0 

0 

0 

24 

76 

200.0 

218473 

.3794 

874 

874 

1391 

1.59 

73 

0 

0 

3 

24 

300.0 

553447 

.9612 

960 

1036 

2370 

2.29 

100 

0 

0 

0 

0 

400.0 

365041 

.6340 

969 

1199 

2789 

2.33 

100 

0 

0 

0 

0 

500.0 

199534 

.3465 

971 

1401 

2651 

1.89 

96 

4 

0 

0 

0 

600.0 

117433 

.2040 

971 

1620 

2511 

1.55 

88 

10 

1 

0 

0 

700.0 

81118 

.1409 

971 

1840 

2598 

1.41 

80 

18 

2 

0 

0 

800.0 

64534 

.1121 

971 

2059 

2728 

1.32 

69 

28 

3 

0 

0 

900.0 

56538 

.  0982 

971 

2279 

2858 

1.25 

59 

37 

4 

0 

0 

1000.0 

52521 

.0912 

971 

2498 

2989 

1.20 

50 

45 

5 

0 

0 

-URSI  - 

-BOGul— 

- 1 

- 1- 

— 

— 

— 

—I 

LAT/LON=  42.5/-71.3  H=1000.0  RZ12=100.0  MMDD:1015  LT:  7.2  SZA=  80.2 

MLA/MLO=  53.8/358.2  DIP=  71.4  F10.7=145.4  DDD:  288  UT:12.0  SDE=  -9.1 


2.2  The  Ionospheric  Conununications  Analysis  And  Prediction  Program  (lONCAP) 

lONCAP  [Lloyd,  et  al. ,  1982]  is  an  improved,  more  versatile  and  flexible  version  of  the  ITS-78 
model  [Barghausen,  et  al. ,  1978].  lONCAP  is  most  useful  for  wave  propagation  using  operating 
frequencies  that  may  be  reflected  by  the  E,  Es,  FI ,  and  F2  layers.  It  is  designed  mainly  for  quiet 
ionospheric  conditions.  lONCAP  does  not  include  the  mid-latitude  trough,  which  exhibits  large 
horizontal  gradients  in  electron  density,  and  does  not  take  into  account  the  effects  of  particle 
precipitation  in  the  auroral  region.  The  model  is  good  for  ±20°  to  ±60°,  but  is  a  poor  predictor 
for  the  equatorial  region  and  the  high  latitude  region. 

The  ITS-78  model  was  developed  by  the  Institute  of  Telecommunication  Sciences,  ESSA, 
Boulder,  Colorado,  primarily  for  the  purpose  of  predicting  long  term  performance  of 
communication  systems  in  the  2  to  30  MHz  frequency  range.  The  important  features  of  ITS-78 
are  the  parameters  for  the  D,  E,  Es,  and  F2  layers  of  the  ionosphere.  lONCAP  contains  several 
significant  improvements  over  ITS-78  [Dandekar,  1982]; 

•  The  description  of  the  ionosphere  is  more  complete. 

•  The  loss  equations  were  supplemented. 

•  Ray  path  geometry  calculations  were  revised. 

•  Loss  statistics  were  revised  to  include  the  effects  of  the  sporadic  E  layer 
and  of  over-the-MUF  (maximum  usable  frequency)  modes. 

•  A  separate  long  path  model  was  developed 

•  The  antenna  gain  package  was  revised. 
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lONCAP  has  thirty  output  options,  which  may  be  divided  into  four  categories:  ionospheric 
description,  antenna  patterns,  MUF  predictions,  and  system  performance  predictions.  Program 
inputs  are  the  date,  universal  time,  geographic  locations  of  the  transmitter  and  receiver,  and 
sunspot  number.  For  system  performance,  some  additional  inputs  are  required,  such  as  radiation 
power  of  the  transmitter,  and  the  signal-to-noise  ratio  of  the  receiver.  External  ionospheric 
parameters  may  also  be  given. 

Table  2  shows  some  sample  lONCAP  output  for  points  along  a  great  circle  path  from  Denver  to 
St.  Louis.  A  full  lONCAP  test  run  required  30  seconds  on  a  66  MHz  486  PC. 


Table  2.  Sample  lONCAP  Output 

JAN  1970  SSN  =  100. 

BOULDER, COLORADO  TO  ST.  LOUIS, MO.  AZIMUTHS  N.  MI.  KM 

40.03  N  105.30  W  -  38.67  N  90.25  W  91.84  281.42  702.6  1301.1 


YE  = 

20.0 

HE  = 

110.0 

HS 

=  110 

.0 

LAT 

LONG 

LMT 

UT 

E 

FI 

Yl 

HI 

FH/2 

F2Z  Y2 

H2 

ES 

MED 

HI 

M3000 

HPF2 

RAT 

ZEN 

2MAX 

MAGL 

39. 6N 

97. 7W 

1.5 

8.0 

.58 

.0 

.0 

.0 

.7 

4.4  70.3 

325.4 

1.3 

1.8 

4.0 

2.97 

325.7 

4.6 

154.7 

66.1 

49. 4N 

39. 6N 

97. 7W 

3.5 

10.0 

.61 

.0 

.0 

.0 

.7 

4.1  75.6 

338.9 

1.1 

1.7 

3.5 

2.89 

339.3 

4.5 

132.9 

66.1 

49. 4N 

39. 6N 

97. 7W 

5.5 

12.0 

1.02 

.0 

.0 

.0 

.7 

3.4  82.2 

335.8 

1.2 

1.7 

3.5 

2.90 

337.8 

4.1 

109.9 

66.1 

49. 4N 

39. 6N 

97. 7W 

7.5 

14.0 

1.95 

.0 

.0 

.0 

.7 

6.0  85.4 

298.3 

1.8 

2.1 

3.9 

3.13 

300.5 

3.5 

88.7 

66.1 

49. 4N 

39. 6N 

97. 7W 

9.5 

16.0 

2.93 

.0 

.0 

.0 

.7 

9.3  93.7 

279.2 

2.6 

2.9 

4.4 

3.26 

281.3 

3.0 

71.7 

66.1 

49. 4N 

39. 6N 

97. 7W 

11.5 

18.0 

3.44 

4.5 

51.4  205.5 

.7 

10.9  105.7 

287.9 

3.1 

3.4 

4.4 

3.16 

295.0 

2.7 

63.0 

66.1 

49. 4N 

39. 6N 

97 .7W 

13.5 

20.0 

3.28 

.0 

.0 

.0 

.7 

11.5  99,5 

304.2 

2.9 

3.2 

4.0 

3.09 

305.9 

3-1 

66.1 

66.1 

49. 4N 

39. 6N 

97. 7W 

15.5 

22.0 

2.43 

.0 

.0 

.0 

.7 

10.8  91.0 

301.4 

2.2 

2.5 

3.5 

3.11 

302.4 

3.3 

79.6 

66.1 

49. 4N 

39. 6N 

97. 7W 

17.5 

.0 

1.38 

.0 

.0 

.0 

.7 

8.7  81.6 

299.2 

1.5 

1.9 

3.5 

3.13 

299.7 

3.7 

99.1 

66.1 

49. 4N 

39. 6N 

97. 7W 

19.5 

2.0 

.74 

.0 

.0 

.0 

.7 

6.0  74.7 

304.5 

1.2 

1.7 

3.9 

3.10 

304.8 

4.1 

121.5 

66.1 

49. 4N 

39. 6N 

97 .7W 

21.5 

4.0 

.55 

.0 

.0 

.0 

.7 

4.3  71.6 

318.7 

1.3 

1.7 

4.4 

3.01 

319.0 

4.5 

144.4 

66.1 

49 .4N 

39. 6N 

97. 7W 

23.5 

6.0 

.56 

.0 

.0 

.0 

.7 

3.8  70.2 

325.6 

1.4 

1.8 

4.4 

2.97 

326.0 

4.6 

162.3 

66.1 

49. 4N 
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2.3  Fully  Analytic  Ionospheric  Model  (FAIM) 

The  Fully  Analytic  Ionospheric  Model  (FAIM)  [Anderson,  et  al ,  1989]  was  designed  to  improve 
over  lONCAP  in  the  equatorial  region.  FAIM  was  created  by  merging  a  modified  version  of  the 
Ching-Chiu  model  [Chiu,  1975]  with  a  semiempirical  low-latitude  ionospheric  model  (SLIM) 
[Anderson,  etal.,  1987]. 

Some  of  the  FAIM  features,  such  as  the  post-sunset  rise  in  the  F  layer  peak  height,  and  the 
equatorial  anomaly  maxima  m  plasma  density  near  +15°  geomagnetic  latitude,  were  designed  to 
approximate  those  in  SLIM.  The  implementation  of  these  features  provides  significant 
improvements  over  previous  comparable  analytic  models. 

The  FAIM  model  will  generate  profiles  of  NO"^,  0"^,  and  total  ion  density  at  one  hour 

increments  of  local  time.  These  profiles  may  be  easily  interpolated  to  any  arbitrary  universal 
time.  As  distributed,  FAIM  provides  the  profiles  over  an  altitude  range  of  90  km  to  1000  km, 
at  10  km  intervals,  although  densities  may  also  be  calculated  for  any  arbitrary  altitude  within  this 
range.  FAIM  is  valid  only  for  geomagnetic  latitudes  between  -40°  and  40°.  Seasonal  variation 
is  controlled  by  specifying  a  month  for  the  run;  no  day  number  may  be  given.  FAIM  has  been 
superseded  by  PIM,  mainly  due  to  this  latitude  restriction. 

The  sample  FAIM  output  shown  in  Table  3  was  generated  for  October,  0  UT,  for  a  flux  of  150. 
The  ion  densities  have  been  converted  to  electron  frequency  using  the  relation  FREQ  =  SQRT 
(DEN/ 12. 4).  Electron  frequencies  are  shown  for  eight  locations  along  a  great  circle  path 
originating  from  Bangor,  Maine.  A  more  extensive  test  run  on  a  40  MHZ  Sparc  IPX  workstation 
using  92  altitude  points,  24  local  times,  and  17  latitude/longitude  pairs  required  79  seconds. 


Table  3.  Sample  FAIM  Output 


Lat 

35.8 

26.8 

Long 

292.2 

292.2 

Range 

1000 

2000 

1000. 

0.402 

0.492 

900. 

0.696 

0.852 

800. 

1.136 

1.391 

700. 

1.796 

2.236 

600. 

2.797 

3.570 

500. 

4.183 

5.448 

400. 

5.658 

7.465 

300. 

5.896 

7.389 

200. 

1.656 

1.606 

100. 

0.635 

0.635 

17.8 

8.8 

-0.2 

292.2 

292.2 

292.2 

3000 

4000 

5000 

0.852 

1.529 

3.626 

1.448 

2.524 

5.529 

2.393 

4.145 

8.309 

3.863 

6.684 

12.045 

6.150 

10.503 

16.072 

9.333 

15.303 

17.579 

12.527 

18.217 

10.962 

11.238 

10.140 

1.555 

1.727 

0.751 

0.635 

0.696 

0.696 

0.696 

-9.2 

-18.2 

-27.2 

292.2 

292.2 

292.2 

6000 

7000 

8000 

4.836 

4.902 

2.965 

6.611 

6.927 

4.701 

8.785 

9.559 

7.432 

11.050 

12.508 

11.486 

12.440 

14.660 

16.841 

10.821 

13.314 

21.540 

5.016 

6.312 

18.766 

0.635 

0.751 

4.803 

0.635 

0.635 

0.751 

0.696 

0.696 

0.751 
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2.4  Parameterized  Ionospheric  Model  (PEM) 

The  Parameterized  Ionospheric  Model  (PIM)  [Whartenby,  1993]  is  the  base  ionospheric  model  on 
which  the  Parameterized  Real-time  Ionospheric  Specification  Model  (PRISM)  [Daniell,  et  al. , 
1995]  operates.  PRISM  uses  data  from  ground-based  and  satellite-based  sources  to  adjust  the 
parameterized  model,  giving  a  near  real-time  specification  of  the  ionosphere.  PIM  is  in  the  public 
domain  and  is  available  for  distribution  to  the  ionospheric  community. 

PIM  is  based  on  four  physical  models: 

•  A  low-latitude  F  layer  model  (LOWLAT)  [Anderson,  1973]. 

•  MIDLAT,  a  mid-latitude  version  of  LOWLAT. 

•  An  E-region  local  chemistry  code  (ECSD)  incorporating  photoelectrons 
using  the  continuous  slowing  down  method  [Jasperse,  1982]. 

•  The  Time  Dependent  Ionospheric  Model  (TDIM)  of  Utah  State  University 
(USU),  a  high  latitude  E  and  F  layer  model  [Schunk,  1988]. 

All  of  these  models  use  the  MSIS-86  thermospheric  model  [Hedin,  1987],  and  are  based  on  a 
tilted  dipole  representation  of  the  geomagnetic  field  and  a  corresponding  magnetic  coordinate 
system.  PIM  is  an  improvement  over  lONCAP  because  it  provides  much  better  results  in  the 
equatorial  region,  mainly  due  to  the  inclusion  of  a  neutral  wind  model.  In  addition,  PIM 
considers  magnetic  activity,  rather  than  assuming  average  conditions. 

PIM  may  produce  either  regional  or  global  ouqjut;  the  latitude/longitude  grid  is  user  selectable. 
The  model  is  valid  over  an  altitude  range  of  90  to  1000  km.  The  following  types  of  ouq)Ut  may 
be  generated: 


•  Electron  density  profiles 

•  Critical  frequencies  and  heights  of  the  E  and  F2  layers 

•  Critical  frequencies  and  heights,  plus  topside  density  and  scale  height,  and 
bottomside  densities. 

The  sample  PIM  output  shown  in  Table  4  was  generated  for  1991,  day  354,  21  UT,  with  a  fixed 
latitude  of  40°.  The  Sun  Spot  Number  was  given  as  130,  and  Kp  was  set  at  2.  The  left  coliunn 
indicates  altitude,  in  kilometers;  east  longitude  is  shown  across  the  top.  A  test  run  on  a  40  MHz 
Sparc  IPX  workstation  calculating  electron  density  profiles  for  92  altitude  points  and  120  latitude 
/longitude  pairs  required  106  seconds. 
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Table  4.  Sample  PIM  Output 


0° 

60° 

120° 

180° 

240° 

300° 

100 

3.89E+03 

8.95E+02 

2.88E+02 

2.68E+03 

1.31E+05 

1.16E+03 

200 

2.06E+03 

2.59E+03 

2.83E+03 

5.18E+05 

5.05E+05 

1. 95E+04 

300 

3.08E+05 

2.22E+05 

1.81E+05 

9.65E+05 

l.llE+06 

7.79E+05 

400 

4.16E+05 

5. 60E+05 

3.13E+05 

4.12E+05 

4.97E+05 

5.60E+05 

500 

2.21E+05 

2. 95E+05 

1.69E+05 

1.98E+05 

2.41E+05 

3.04E+05 

600 

1.06E+05 

1.28E+05 

8.02E+04 

1.12E+05 

1.35E+05 

1.70E+05 

700 

5.47E+04 

5. 69E+04 

4.00E+04 

7.46E+04 

8.98E+04 

l.llE+05 

800 

2.92E+04 

2. 62E+04 

2.06E+04 

5.20E+04 

6.22E+04 

7.48E+04 

900 

1.60E+04 

1.23E+04 

1.08E+04 

3.73E+04 

4.44E+04 

5.21E+04 

1000 

8.97E+03 

5. 95E+03 

5.81E+03 

2.72E+04 

3.22E+04 

3.70E+04 

2.5  General  Electric  Over-the-Horizon  (GE/OTH)  Ionospheric  Model 


The  General  Electric  Over-the-Horizon  (GE/OTH)  model  is  a  revision  of  the  AFGWC  polar 
model  [AFGWC,  1982],  developed  by  the  U.S.  Air  Force  Global  Weather  Central.  This  model, 
also  known  as  the  AN/FPS  118  -  Analytical  Model  Specification  [Millman,  et  al,  1988],  was 
developed  for  frequency  management  of  the  Over-the-Horizon  Backscatter  (OTH-B)  radar,  and 
for  target  coordinate  registration  and  propagation  assessments. 

The  GE/OTH  model  uses  the  ITS-78  ionospheric  model  [Barghausen,  1978]  as  a  starting  point. 
The  basic  model  is  driven  by  an  effective  sunspot  number  computed  from  foFj  data  for  the  past 
5  days  from  a  network  of  50  ionospheric  stations.  Using  Kp  and  the  effective  auroral  oval 
parameter  Qg,  the  auroral  E  and  F  layer  enhancements  and  the  mid-latitude  F  layer  trough  electron 
density  depletion  are  superposed  on  the  ITS-78  model.  Ionospheric  parameters  obtained  by  the 
automated  Air  Weather  Service  (AWS)  Digital  Ionospheric  Sounding  System  (DISS)  network  are 
used  for  real-time  updating  of  the  model  at  the  radar  sites. 

The  software  for  the  GE/OTH  model  is  somewhat  more  complicated  than  that  of  other  ionospheric 
models.  Presently,  the  program  is  available  only  for  the  VAX/VMS  operating  system.  Two 
separate  connections  are  required  -  a  "master"  program,  which  runs  the  actual  model,  and  a 
"driver",  which  provides  the  user  interface.  The  two  processes  are  linked  through  interprocess 
communication  and  a  shared  memory  area.  A  full  program  run  takes  approximately  90  seconds 
on  a  VAX  7000. 

The  GE/OTH  model  requires  a  "grid"  file  containing  a  series  of  latitude/longitude  pairs.  The  user 
must  provide  Kp,  Q,  sunspot  number,  date,  and  time;  these  values  are  entered  through  the  driver 
program.  Table  5  presents  some  typical  output  from  the  GE/OTH  model.  The  ground  points 
were  traced  along  a  170°  azimuth  from  the  East  Coast  Radar  System  (ECRS)  in  Bangor,  Maine. 
The  run  was  made  for  December  20,  1992,  using  Kp  and  Q  of  2,  and  a  sunspot  number  of  75. 
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Table  5.  Sample  GE/OTH  Model  Output 


Range 

GGLAT 

GGLON 

GMLAT 

GMLON 

foF2 

HmF2 

YinF2 

foFl 

HmFl 

YitiFl 

foE 

HrnE 

YmE 

16 

44.8 

292.0 

55.9 

10.0 

4.10 

317.9 

69.3 

0.00 

0.0 

0.0 

0.60 

120.0 

21.8 

998 

35.9 

293.9 

46.9 

11.2 

4.27 

318.4 

69.1 

0.00 

0.0 

0.0 

0.60 

120.0 

21.8 

1999 

27.0 

295.4 

38.0 

11.7 

4.35 

318.0 

75.3 

0.00 

0.0 

0.0 

0.60 

120.0 

21.8 

3000 

18.1 

296.6 

29.3 

11.8 

4.64 

315.7 

79.6 

0.00 

0.0 

0.0 

0.60 

120.0 

21.8 

4001 

9.1 

297.8 

20.9 

11.6 

6.35 

312.1 

80.9 

0.00 

0.0 

0.0 

0.60 

120.0 

21.8 

5001 

0.2 

298.9 

12.6 

11.3 

10.15 

314.1 

83.1 

0.00 

0.0 

0.0 

0.60 

120.0 

21.8 

6002 

-8.7 

300.0 

4.3 

10.9 

9.69 

345.7 

91.5 

0.00 

0.0 

0.0 

0.60 

120.0 

21.8 

7003 

-17.6 

301.1 

-5.1 

10.6 

9.06 

388.2 

102.9 

0.00 

0.0 

0.0 

0.60 

120.0 

21.8 

8004 

-26.6 

302.4 

-13.9 

10.6 

10.02 

385.2 

103.3 

0.00 

0.0 

0.0 

0.60 

120.0 

21.8 

9005 

-35.5 

303.9 

-22.3 

11.2 

8.99 

375.2 

103.8 

0.00 

0.0 

0.0 

0. 60 

120.0 

21.8 

10005 

-44  .4 

305.7 

-30.5 

12.4 

8.55 

369.2 

105.1 

1.01 

231.5 

57.9 

0.60 

120.0 

21.8 

11006 

-53.2 

308.2 

-38.7 

14.5 

8.20 

350.0 

98.5 

2.49 

227.4 

56.9 

1.83 

120.0 

21.8 

12007 

-62.0 

312.1 

-47.0 

17.5 

7.50 

330.5 

92.4 

3.65 

223.3 

55.8 

2.34 

120.0 

21.8 

13008 

-70.6 

319.0 

-55.4 

21.8 

6.65 

310.2 

90.1 

4.00 

219.2 

54.8 

2.64 

120.0 

21.8 

2.6  The  RAYLAB  Ray  Tracing  Program 

RAYLAB  is  a  two-dimensional  ray  tracing  program  developed  by  Chris  Coleman  of  the  Australian 
Defence  Science  and  Technology  Organization  (DSTO).  RAYLAB  was  developed  concurrently  with 
HASEL  [Coleman,  1993],  a  three-dimensional  ray  tracing  program. 


A  two-dimensional  ray  trace  assumes  that  the  rays  are  confined  to  a  fixed  azimuthal  path  emanating 
from  the  radar.  Three-dimensional  models  more  closely  approximate  actual  ray  behavior;  however, 
the  “drift:”  in  and  out  of  the  azimuthal  plane  is  generally  small  compared  to  the  distance  traveled  from 
the  radar.  The  two-dimensional  model  may  be  constructed  with  a  substantially  reduced  set  of 
equations,  resulting  in  a  significant  increase  in  program  speed. 

RAYLAB  is  designed  to  ray  trace  in  an  ionosphere  with  a  spatially  varying  electron  density.  A 
Runge-Kutta-Fehlberg  scheme  is  used  to  solve  a  simplified  version  of  the  Haselgrove  form  of  the  ray 
equations.  In  essence,  this  procedure  is  a  simplified  version  of  the  technique  used  in  HASEL.  Two 
slightly  different  versions  of  RAYLAB  have  been  distributed.  One  uses  a  simple  Chapman  layer 
ionospheric  model,  while  the  other  incorporates  the  FAIM  model.  In  addition,  a  third  version  has 
been  adapted  to  use  ionospheric  predictions  from  PIM. 

A  portion  of  the  output  from  a  typical  RAYLAB  run  is  shown  in  Table  6.  The  ionosphere  was 
generated  with  PIM,  for  1991,  day  354, 1  UT,  with  a  Sunspot  Number  of  109.4  and  Kp  of  2.  Rays 
were  traced  along  a  1 10°  azimuth  from  Bangor,  Maine  (44.8°  N,  -68°  E),  at  a  frequency  of  10  MHz. 
The  data  given  in  the  table  show  altitude  and  range  for  a  ray  with  a  starting  elevation  of  4°,  out  to  a 
range  of  6000  km.  The  complete  run,  calculating  rays  for  16  different  starting  elevations,  required 
1.06  seconds  on  a  Silicon  Graphics  workstation  equipped  with  a  150  MHZ  R4400  processor. 
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Table  6.  Sample  RAYLAB  Output 


Alt 

Range 

Alt 

Range 

0.0 

0.0 

248.8 

1538.7 

3.7 

49.8 

242.2 

1585.0 

7.8 

99.6 

232.6 

1631.5 

12.2 

149.4 

221.0 

1678.1 

17.1 

199.0 

208.8 

1724.9 

22.3 

248.6 

196.6 

1771.9 

27.9 

298.1 

184.8 

1819.0 

33.9 

347.5 

173.4 

1866.3 

40.3 

396.8 

162.3 

1913.8 

47.1 

446.0 

151.7 

1961.4 

54.2 

495.1 

141.4 

2009.2 

61.8 

544.1 

131.5 

2057.1 

69.7 

592.9 

122.0 

2105.2 

78.0 

641.7 

112.7 

2153.4 

86.7 

690.3 

103.6 

2201.8 

89.9 

708.0 

95.0 

2250.3 

95.5 

737.6 

86.7 

2298.9 

104.8 

785.9 

78.8 

2347.6 

114.5 

834.1 

71.3 

2396.5 

124.6 

882.2 

64.2 

2445.5 

135.0 

930,1 

57.4 

2494.5 

145.8 

977.8 

51.0 

2543.7 

156.9 

1025.4 

45.1 

2593.0 

168.3 

1072.8 

39.4 

2642.3 

180.0 

1120.1 

34.2 

2691.8 

192.2 

1167.1 

29.4 

2741.3 

204.7 

1214.0 

24.9 

2790.9 

217.3 

1260.8 

20.9 

2840.6 

229.3 

1307.3 

17.2 

2890.3 

239.6 

1353.7 

13.9 

2940.0 

247.2 

1400.0 

11,0 

2989.9 

251.3 

1446.3 

8.5 

3039.7 

251.9 

1492.5 

Alt 

Range 

Alt 

Range 

6.4 

3089.6 

134.8 

4626.8 

4.7 

3139.5 

145.0 

4674.6 

3.3 

3189.5 

155.5 

4722.4 

2.4 

3239.5 

166.3 

4769.9 

1.9 

3289.4 

177.4 

4817.3 

1.7 

3339.4 

188.8 

4864.6 

1.9 

3389.4 

197.0 

4898.4 

2.6 

3439.4 

204.1 

4928.4 

3.6 

3489.4 

209.2 

4952.1 

5.0 

3539.3 

214.0 

4979.0 

6.8 

3589.2 

217.3 

5005.6 

9.0 

3639.1 

218.7 

5035.3 

11.6 

3689.0 

216.2 

5082.5 

14.5 

3738.8 

209.4 

5129.8 

17.9 

3788.5 

200.8 

5177.3 

21.7 

3838.2 

192.0 

5224.9 

25.8 

3887.9 

183.4 

5272.7 

30.3 

3937.5 

175.2 

5320.6 

35.2 

3987.0 

167.4 

5368.7 

40.5 

4036.4 

160.0 

5416.8 

46.2 

4085.7 

152.9 

5465.1 

52.3 

4135.0 

146.2 

5513.4 

58.7 

4184.1 

139.9 

5561.9 

65.6 

4233.2 

133.7 

5610.4 

72.8 

4282.1 

127.8 

5659.1 

80.4 

4331.0 

122.3 

5707.8 

88.3 

4379.7 

117.1 

5756.6 

90.8 

4394.1 

112.4 

5805.4 

99.3 

4442.6 

108.2 

5854.4 

106.4 

4482.3 

104.5 

5903.4 

115.5 

4530.6 

101.0 

5952.4 

125.0 

4578.8 

97.6 

6001.5 

2.7  The  HASEL  Ray  Tracing  Program 


BASEL  [Coleman,  1993]  is  a  three-dimensional  ray  tracing  program  developed  by  Chris  Coleman 
of  the  Australian  Defence  Science  and  Technology  Organization  (DSTO).  HASEL  was  developed 
concurrently  with  RAYLAB,  a  two-dimensional  ray  tracing  program. 


Unlike  a  two-dimensional  ray  trace,  which  assumes  that  the  rays  are  confined  to  a  fixed  azimuthal 
path,  three-dimensional  models  closely  approximate  actual  ray  behavior.  Consequently,  the  three- 
dimensional  model  is  constructed  with  a  much  more  complex  set  of  equations,  resulting  in  a 
significant  increase  in  program  run-time. 

The  procedure  is  based  on  a  numerical  solution  to  the  Haselgrove  ray  tracing  equations.  These 
are  ordinary  differential  equations,  which  are  solved  using  a  Runge-Kutta-Fehlberg  scheme.  The 
Earth’s  magnetic  field  is  represented  by  a  tilted  dipole.  There  are  three  options  for  representing 
the  ionosphere; 
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•  The  ionosphere  can  be  defined  using  an  internal  function  subroutine,  ELDEN.  This  is  a  simple 
analytic  ionospheric  model. 


•  A  description  of  the  ionosphere  may  be  given  in  terms  of  layer  parameters,  defined  on  a  regular 
(that  is,  even  increments)  geographic  latitude/longitude  grid.  The  three  layers  (E,  F],  and  Fj)  are 
defined  by  their  heights,  plasma  frequencies  at  these  heights,  and  layer  thicknesses.  The  electron 
density  profile  is  then  composed  of  Chapman  layers. 

•  The  third  option  is  similar  to  the  previous  method,  but  the  Chapman  layers  are  replaced  with  a 
set  of  samples  above  each  geographical  grid  point.  As  with  the  geographic  sampling,  the  altitude 
sampling  must  also  be  uniform;  a  three-dimensional  interpolation  scheme  is  then  used  to  determine 
electron  densities.  This  is  the  most  efficient  of  the  three  ionospheric  representations. 


A  portion  of  the  output  from  a  typical  BASEL  run  is  shown  in  Table  7.  This  ray  trace  was  made 
using  the  ELDEN  function  to  represent  the  ionosphere,  with  the  magnetic  field  effects  active,  and 
medium  tolerance.  Rays  were  evaluated  for  all  combinations  of  four  bearings,  six  initial  elevation 
angles,  and  three  frequencies.  The  complete  run  required  66  seconds  on  a  Silicon  Graphics 
workstation  equipped  with  a  150  MHZ  R4400  processor. 


Table  7.  Sample  BASEL  Output 


Doppler 

Shift 

Ground 

Range 

Group 

Path 

Phase 

Path 

Max 

Height 

Init 

Elev 

Final 

Elev 

Final 

Long 

Final 

Lat 

Bear 

Freq 

Mode 

-0.84 

1882.2 

1975.7 

1956.3 

219.3 

10.0 

10.0 

295.3 

28.1 

170.0 

10.0 

F20 

-1.69 

3769.4 

3956.3 

3917.6 

219.5 

10.0 

9.9 

297.7 

11.3 

170.0 

10.0 

F20F20 

-1.68 

3751.4 

3938.0 

3898.8 

217.2 

10.0 

10.0 

297.6 

11.4 

170.0 

10.0 

F20F2 

-2.53 

5642.5 

5922.7 

5864.2 

217.1 

10.0 

9.9 

299.8 

-5.4 

170.0 

10.0 

F20F20F2X 

-2.53 

5661.0 

5941.6 

5883.6 

219.5 

10.0 

9.9 

299.8 

-5.6 

170.0 

10.0 

F2OF20F20 

-2.52 

5622.7 

5902.6 

5843.5 

217.2 

10.0 

10.0 

299.7 

-5.2 

170.0 

10.0 

F20F2XF2X 

-2.53 

5641.2 

5921.5 

5862.9 

219.6 

10.0 

9.9 

299.8 

-5.4 

170.0 

10.0 

F20F2XF20 

-0.84 

1865.7 

1958.8 

1939.0 

217.3 

10.0 

10.0 

295.3 

28.2 

170.0 

10.0 

F2X 

-1.69 

3749.7 

3936.5 

3897.3 

219.5 

10.0 

10.0 

297.6 

11.5 

170.0 

10.0 

F2XF20 

-1.68 

3731.8 

3918.2 

3878.4 

217.2 

10.0 

10.0 

297.6 

11.6 

170.0 

10.0 

F2XF2X 

-2.53 

5619.8 

5899.8 

5840.8 

217.2 

10.0 

10.0 

299.7 

-5.2 

170.0 

10.0 

F2XF20F2X 

-2.53 

5638.2 

5918.6 

5860.1 

219.5 

10.0 

10.0 

299.8 

-5.4 

170.0 

10.0 

F2XF20F20 

-2.52 

5600.0 

5879.7 

5820.1 

217.2 

10.0 

10.0 

299.7 

-5.0 

170.0 

10.0 

F2XF2XF2X 

-2.53 

5618.5 

5898.6 

5839.5 

219.6 

10.0 

10.0 

299.7 

-5.2 

170.0 

10.0 

F2XF2XF20 

3.  RADIATION  BELT  FLUX  AND  DOSE  MODELS 

The  radiation  belts  are  regions  within  the  earth’s  magnetosphere  where  energetic  ions  and 
electrons  are  trapped  by  the  earth’s  magnetic  field  at  radial  distances  between  1.1  and  10  earth 
radii  [Spjeldvik  and  Rothwell,  1985].  The  energy  of  these  charged  particles  range  between  1  and 
200  MeV  and  are  sufficient  to  penetrate  inadequately  shielded  spacecraft,  where  they  can  cause 
serious  damage  to  electronic  components  operating  on  spacecraft,  including  electronic  discharges, 
degradation,  and  alteration  of  memory  (single  event  upsets). 
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The  flux  and  dose  models  described  here  are  based  upon  measurements  from  the  Combined 
Release  and  Radiation  Effects  Satellite  (CRRES)  mission,  flown  in  1991-1992,  and  from 
measurements  from  various  satellites  (the  NASA  models)  during  the  period  1958-1970. 

The  CRRES  models  are  incorporated  in  the  CRRESRAD,  CRRESPRO  and  CRRESELE  software 
packages.  CRRESRAD  contains  the  CRRES  dose  model,  CRRESPRO  the  proton  flux  model,  and 
CRRESELE,  the  electron  flux  model. 

The  Combined  Release  and  Radiation  Effects  Satellite  (CRRES)  provided  the  first  comprehensive 
study  of  the  radiation  belts  since  the  NASA  models.  CRRES  was  launched  on  25  July  1990  and 
provided  measurements  in  a  geosynchronous  transfer  orbit  (perigee  altitude  350  km,  apogee 
altitude  33000  km,  inclination  18°,  orbital  period  9.85  hr)  until  12  October  1991.  CRRES 
contained  a  wide  array  of  instruments  to  measure  radiation  belt  fluxes  and  dosage  and  their  effects 
on  operating  components  of  satellites. 

The  CRRES  models  software  packages  are  available  from  PL  in  the  form  of  PC  executable  codes 
and  related  binary  data  files.  The  programs  include  an  easy-to-use  graphics  visualization  package 
to  facilitate  examination  of  the  data  by  the  user  and  to  compute  relevant  physical  properties  of  the 
radiation  at  given  locations  in  space  or,  where  relevant,  integrated  over  satellite  orbits.  A  custom 
version  of  the  LOKANGL  orbital  prediction  program  is  provided  as  part  of  the  CRRES  model 
packages  for  computation  of  flux  or  dose  over  a  specified  orbit.  The  documentation  provides  a 
description  of  the  model  binary  data  files  that  is  sufficient  to  facilitate  the  use  of  the  models  in 
custom  user  written  applications.  The  sample  output  data  provided  in  this  section  is  based  upon 
custom  applications  written  by  Radex  based  upon  the  binary  data  model  files.  Sample  output  is 
provided  here  only  for  the  CRRESRAD  model. 

The  NASA  flux  models  are  AP8MAX,  AP8MIN  (for  protons),  AE8MAX,  and  AE8MIN  (for 
electrons)  where  MAX  and  MIN  refer  to  solar  maximum  and  solar  minimum  respectively.  Model 
data  files  and  FORTRAN  source  code  are  provided  for  flux  and  dosage  computations. 


3.1  CRRESRAD,  The  CRRES  Dose  Model 

Radiation  dosage  is  a  measure  of  the  energy  deposited  by  radiation  in  a  material  per  unit  mass  of 
the  material.  It  is  usually  given  in  rads,  where  1  rad  =  100  ergs  per  gram.  On  CRRES,  the  dose 
was  measured  by  the  Space  Radiation  Dosimeter  [Hardy,  et  al.,  1985;  Morel,  et  al.,  1989], 
which  consisted  of  4  silicon  detectors  shielded  by  hemispheres  of  aluminum  with  thicknesses  of 
82.5,  232.5,  457.4,  and  886.5  mils  (1  mil  =  1/1000  inch).  These  correspond  to  energy 
thresholds  of  20,  35,  50,  and  75  MeV  for  protons,  and  1,  2.5,  5,  and  10  MeV  for  electrons.  For 
each  detector,  dose  was  accumulated  in  two  channels,  LOLET  and  HILET,  corresponding  to 
single  event  energy  depositions  of  0.01  -  1  MeV  and  1-10  MeV,  respectively.  HILET  dose  is 
caused  principally  by  protons  of  energies  20-100  MeV,  and,  in  major  magnetic  storms,  electrons 
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above  5  MeV.  LOLET  dose  is  caused  by  electrons,  bremsstrahlung,  and  protons  above  100  MeV. 

A  major  geomagnetic  storm  in  late  March,  1991,  so  profoundly  changed  the  magnetosphere 
[Mullen,  et  al. ,  1991;  Blake,  et  al. ,  1992]  that  two  dose  models  were  developed  for  each  shield 
thickness  and  dose  channel;  one  for  quiet  conditions  ("Quiet"  model),  based  on  data  before  the 
storm,  and  the  other  for  active  conditions  ("Active"  model,  based  on  data  after  the  storm 
[Gussenhoven,  et  al,  1992;  Kems  and  Gussenhoven,  1992].  The  total  number  of  models 
available  is  36:  for  each  of  the  4  thicknesses  there  are  models  for  each  channel  (HILET  and 
LOLET)  and  a  model  for  combined  HILET  and  LOLET.  Each  of  these  is  given  for  quiet  and 
active  conditions,  and  the  entire  CRRES  mission.  Each  model  specifies  dose  rate  in  bins  of  L  and 
B/Bq.  L  is  the  Mcllwain  shell  parameter  [Mcllwain,  1961]  which  for  a  dipole  field  is  the 
equatorial  crossing  distance  of  the  magnetic  field  line  through  the  point  of  interest.  B  is  the 
magnetic  field  intensity  at  the  point  of  interest,  while  Bq  is  the  intensity  at  the  equatorial  crossing 
of  the  field  line.  The  bin  boundaries  are  defined  by: 

Lj  =  0.95  +  0.05  I  Re,  1  <  =  I  <  =  141 

B/Bo  =  1.000,  1.004,  1.020,  1.046,  1.085,  1.140,  1.200,  1.300,  1.400,  1.520,  1.690,  1.880, 
2.100,  2.400,  2.730,  3.130,  3.670,  4.350,  5.020,  6.100,  7.410 

These  models  have  been  incorporated  into  the  CRRESRAD  model  [Kems  and  Gussenhoven, 
1992],  which  predicts  the  radiation  dose  received  in  a  specified  orbit. 

The  test  case  shown  in  Table  8  provides  sample  dose  rates  (rads  per  sec)  at  the  magnetic  equator 
(B/Bq  =  1)  for  the  various  channels,  conditions,  thicknesses  and  L  values  indicated.  A  test  on  a 
40  MHZ  Sparc  workstation  making  10000  calls  to  the  dose  model  routine  required  6.5  seconds. 
The  program  which  was  used  to  generate  Table  8  was  written  at  Radex  based  upon  information 
about  the  binary  data  files  provided  in  the  program  documentation  [Kems  and  Gussenhoven, 
1992]. 
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Table  8.  Sample  CRRES  Dose  Model  Output 


Channel 

Combined 

Combined 

Combined 

Combined 

Combined 

Combined 

Combined 

Combined 

Combined 

Combined 

Combined 

Combined 

Hilet 

Hilet 

Hilet 

Hilet 

Hilet 

Hilet 

Hilet 

Hilet 

Hilet 

Hilet 

Hilet 

Hilet 

Lolet 

Lolet 

Lolet 

Lolet 

Lolet 

Lolet 

Lolet 

Lolet 

Lolet 

Lolet 

Lolet 

Lolet 


Condition 

Combined 

Combined 

Combined 

Combined 

Quiet 

Quiet 

Quiet 

Quiet 

Active 

Active 

Active 

Active 

Combined 

Combined 

Combined 

Combined 

Quiet 

Quiet 

Quiet 

Quiet 

Active 

Active 

Active 

Active 

Combined 

Combined 

Combined 

Combined 

Quiet 

Quiet 

Quiet 

Quiet 

Active 

Active 

Active 

Active 


Thickness 

(Mils) 

82.5 

232.5 

457.5 

886.5 

82.5 

232.5 

457.5 

886.5 

82.5 

232.5 

457.5 

886.5 

82.5 

232.5 

457.5 

886.5 

82.5 

232.5 

457.5 

886.5 

82.5 

232.5 

457.5 

886.5 

82.5 

232.5 

457.5 

886.5 

82.5 

232.5 

457.5 

886.5 

82.5 

232.5 

457.5 

886.5 


1.50 

1.834E-02 
4.048E-03 
2 .171E-03 
1.199E-03 
1.807E-02 
4.031E-03 
2.150E-03 
1.167E-03 
1.869E-02 
4.080E-03 
2.197E-03 
1.233E-03 
1.598E-02 
3 .356E-03 
1.538E-03 
8 .064E-04 
1.582E-02 
3.356E-03 
1.537E-03 
8.003E-04 
1.619E-02 
3 .369E-03 
1.542E-03 
8.131E-04 
2 .360E-03 
6.914E-04 
6.326E-04 
3 .921E-04 
2 .246E-03 
6.748E-04 
6.128E-04 
3.664E-04 
2 .499E-03 
7.108E-04 
6 .555E-04 
4.197E-04 


L  (RE)--- 

2,50 

1.722E-03 
8.409E’-04 
4.565E-04 
1.146E-04 
4.473E-05 
2 .676E-05 
8.410E-06 
4.330E-06 
2.995E-03 
1.469E-03 
7.869E-04 
1.962E-04 
4.532E-04 
1.913E-04 
6.583E-05 
1.042E-05 
2.162E-05 
1.071E-05 
3.470E-06 
2 .710E-06 
7 .818E-04 
3 .312E-04 
1.150E-04 
1.614E-05 
1.269E-03 
6.497E-04 
3 .907E-04 
1.042E-04 
2.311E-05 
1.605E-05 
4.940E-06 
1.620E-06 
2.213E-03 
1.138E-03 
6.719E-04 
1.801E-04 


4.50 

1.401E-02 
1.995E-04 
2.806E-05 
1.289E-05 
4.849E-03 
4.003E-05 
1.013E-05 
5.810E-06 
2.426E-02 
3,900E-04 
4.779E-05 
2.111E-05 
5.060E-06 
3 .290E-06 
O.OOOE+00 
O.OOOE+00 
O.OOOE+00 
l,600E-06 
O.OOOE+00 
O.OOOE+00 
6.540E-06 
5.180E-06 
O.OOOE+00 
4.600E-07 
1.401E~02 
1.962E~04 
2.806E-05 
1.289E-05 
4.849E-03 
3.843E-05 
1.013E-05 
5.810E-06 
2.425E-02 
3.849E-04 
4 .779E-05 
2.065E-05 
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3.2  CRRESPPRO,  The  CRRES  Proton  Model 

The  CRRES  proton  flux  model  CRRESPRO  [Meffert  and  Gussenhoven,  1994\  is  based  upon 
measurements  made  by  the  Proton  Telescope  (PROTEL),  which  was  used  to  measure  proton 
energy  and  flux  on  the  CRRES  satellite  in  the  range  1-100  MeV.  PROTEL  is  a  sophisticated 
proton  detector  designed  to  operate  in  a  hostile  environment,  where  it  was  subjected  to  a  high 
density  flux  of  electrons,  protons  and,  to  a  lesser  extent,  heavier  ions,  particularly  while  passing 
through  the  radiation  (Van  Allen)  belts.  In  addition  to  using  passive  shielding  techniques  and 
magnetic  deflection  of  electrons,  it  used  an  array  of  silicon  particle  detectors  together  with  an  on¬ 
board  processor,  that  used  coincidence  and  anti-coincidence  in  the  detector  array  (detection  logic) 
to  reduce  false  counts.  A  detailed  description  of  the  PROTEL  instrument  is  provided  in  {Lynch, 
etal,  1989]  and  [Violet,  etal.,  1992]. 

PROTEL  consisted  of  two  detector  instruments,  the  low  energy  head  [LEH]  (1-10  MeV)  and  the 
high  energy  head  [HEH]  (6  -  100  MeV)  and  a  dedicated  processor  that  processes  the  raw  data 
from  the  detectors,  and  at  1  second  intervals,  hands  off  the  reduced  data  to  the  satellite’s  telemetry 
system.  The  reduced  data  consist  of  counts  in  24  energy  channels  (8  for  LEH,  16  for  HEH) 
spaced  logarithmically  in  the  1  - 100  MeV  energy  interval,  together  with  environmental  data  and 
raw  counts  for  the  solid  state  particle  detectors.  Flux  data  files  are  provided  for  22  energy 
channels  centered  at  the  following  proton  energies:  1.5,  2.1,  2.5,  2.9,  3.6,  4.3,  5.7,  6.8,  8.5, 
9.7,  10.7,  13.2,  16.9,  19.4,  26.3,  30.9,  36.3,  42.3,  47.5,  57.0,  67.5  and  92.9  MeV.  For 
flux/fluence  calculations  the  CRRESPRO  model  software  omits  the  3.6  and  16.5  MeV  channels 
because  of  overlap. 

As  in  the  case  of  the  CRRES  dose  model,  two  flux  models,  "Quiet"  and  "Active",  are  provided 
for  CRRESPRO,  corresponding  to  data  collected  before  and  after  the  March  1991  storm. 

Each  model  specifies  dose  rate  in  bins  of  L  and  B/Bq.  L  is  the  Mcllwain  shell  parameter 
[Mcllwain,  1961]  which  for  a  dipole  field  is  the  equatorial  crossing  distance  of  the  magnetic  field 
line  through  the  point  of  interest.  B  is  the  magnetic  field  intensity  at  the  point  of  interest,  while 
Bq  is  the  intensity  at  the  equatorial  crossing  of  the  field  line.  The  bin  boundaries  are  defined  by: 

Li  =  0.95  +  0.05  i  Re,  1  <  =  i  <  =  90 

B/Bo  =  1.000,  1.004,  1.020,  1.046,  1.085,  1.140,  1.200,  1.300,  1.400,  1.520,  1.690,  1.880, 
2.100,  2.400,  2.730,  3.130,  3.670,  4.350,  5.020,  6.100,  7.410,  9.088,  11.29,  14.22,  18.16, 
23.56,  31.07,  41.47,  57.23,  80.25,  115.4,  170.7,  260.7,  413.4  and  684.6.  These  values  cover 
approximately  ±68°  magnetic  latitude  in  a  dipole  field. 
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3.3  CRRESELE,  The  CRRES  Electron  Model 

The  CRRES  Electron  model  CRRESELE  [Brautigam  and  Bell,  1995]  is  based  upon  measure¬ 
ments  made  by  the  High  Energy  Electron  Fluxmeter  (HEEF)  instrument  on  the  CRRES  satellite. 
HEEF  measured  the  electron  flux  of  1-10  MeV  in  ten  differential  flux  number  channels.  Eight 
models  are  provided,  six  of  which  correspond  to  average  flux  for  different  levels  of  solar  activity , 
and  two  corresponding  to  observed  daily  average  and  daily  maximum  flux. 

The  L-bin  boundaries  are  defined  by: 

Li  =  2.45  +  0.05i  Rg,  1  <  =  i  <  =  86 

The  B/Bq  bins  are  the  same  as  in  the  CRRESPRO  model. 


3.4  NASA  Radiation  Belt  Trapped  Particle  Flux  Models  AP-8  And  AE-8 

The  empirical  models  AP-8  and  AE-8  describe  radiation  belt  trapped  proton  and  electron  fluxes, 
respectively,  in  the  ion  energy  range  0. 1-400  MeV  and  electron  energy  range  0.04-7  MeV.  They 
are  based  on  data  from  1958  to  1970.  They  are  described  briefly  by  Jordan  [1989]  with  further 
details  given  in  the  references  therein.  The  models  give  integral  or  differential  onmidirectional 
flux  as  a  functions  of  L  [Mdlwain,  1961]  and  B/Bq,  where  B  is  the  magnetic  field  magnitude  at 
the  point  of  interest,  and  Bq  is  the  equatorial  B  value,  given  by: 

B^igauss)  =  0.311653/L^ 


For  each  model  there  is  a  solar  maximum  (AP8MAX,  AE8MAX)  and  solar  minimum  (AP8MIN, 
AE8MIN)  version.  Other  variations  found  in  the  data,  including  those  related  to  magnetic  storms 
and  local  time,  are  averaged  in  these  models.  The  data  supporting  the  earliest  models  was 
contaminated  by  the  Starfish  nuclear  detonation  (9  July  1962),  but  this  has  been  removed  in  the 
present  models. 

Table  1  of  Sawyer  and  Vette  [1976]  lists  the  data  sources  for  AP-8.  These  authors  found  that  in 
the  inner  zone  (L  <  2.5  Rg)  the  fluxes  were  generally  stable  except  for  responses  to  major  storms, 
the  decay  of  the  starfish  population,  and  solar  cycle  dependence  at  altitudes  less  than  1000  km. 
In  the  outer  zone,  nonadiabatic  storm-related  flux  changes  were  observed,  with  recovery  time 
constants  ranging  from  months  to  years.  Modeling  of  the  storm-related  variations  was  not 
attempted  because  of  the  difficulty  in  generalizing  the  details  of  specific  storms.  The  model  solar 
cycle  variation  was  derived  from  the  satellite  Azur  data.  Although  earlier  models  were  generated 
by  analytic  fits  to  power  law  or  exponential  spectra,  the  present  AP-8  models  are  purely  tabular 
in  form,  because  of  the  difficultly  in  fitting  the  entire  0.1-400  MeV  range  to  such  simple  forms. 
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The  AE-8  models  were  derived  from  a  combination  of  the  AE-4,  AE-5  and  AE-6  models  with  new 
data  from  satellites  OVl-3,  OVl-19,  AZUR,  and  ATS-6.  The  AE-4  model  was  derived  from  the 
measurements  of  23  instruments  flown  on  11  satellites  from  1959  to  1968.  The  AE-5  model  was 
derived  from  data  between  December,  1964,  and  December,  1967,  covering  the  inner  zone  only 
(L  =  1 .2  -  2.8  Re).  The  form  used  for  AE-8,  called  AE-5  1975  projected,  models  solar  minimum 
conditions.  AE-6  is  used  to  model  solar  maximum.  Both  AE-5  1975  projected  and  AE-6  have 
had  the  Starfish  contribution  removed.  Unlike  earlier  models,  which  were  tabulated  separately 
for  inner  and  outer  zones,  the  AE-8  model  tables  are  combined,  just  like  the  AP-8  models.  The 
AP-8  and  AE-8  tables  use  common  format,  permitting  them  to  be  used  in  the  same  software. 

The  sample  ouQ)ut  in  Table  9  was  generated  for  1  MeV  ions  (AP8MAX)  and  electrons 
(AE8MAX)  for  solar  maximum  conditions. 


A  test  on  a  40  MHz  Sparc  IPX  workstation  making  10000  calls,  each  reading  in  a  model  table 
from  a  data  file,  required  340  sec.  A  second  test,  also  making  10000  calls,  but  reading  in  a  table 
only  on  the  first  call,  required  0.4  sec. 


Table  9.  Sample  NASA  Model  Output  (1  MeV  Particles) 

AP8MAX  MODEL 


B/BO 


L 

1.0 

o 

CM 

4.0 

6.0 

GO 

O 

10.0 

30.0 

1.0 

O.OOOE+00 

O.OOOE+OO 

O.OOOE+OO 

O.OOOE+OO 

O.OOOE+OO 

O.OOOE+OO 

O.OOOE+OO 

2.0 

6.928E+00 

5.904E+00 

4.935E+00 

3.333E+00 

O.OOOE+OO 

O.OOOE+OO 

O.OOOE+OO 

3.0 

7.580E+00 

6.764E+00 

6.203E+00 

5.884E+00 

5.646E+00 

5.444E+00 

O.OOOE+OO 

4.0 

6.493E+00 

5.835E+00 

5.390E+00 

5.159E+00 

4.985E+00 

4.870E+00 

4.171E+00 

5.0 

4.527E+00 

3.599E+00 

3. 080E+00 

2.848E+00 

2.683E+00 

2.541E+00 

1.977E+00 

6.0 

3.415E+00 

2.491E+00 

2.024E+00 

1.811E+00 

1. 649E+00 

1.543E+00 

9.740E-01 

7.0 

2.669E+00 

1.826E+00 

1.369E+00 

1.174E+00 

1.043E+00 

9.305E-01 

3.999E-01 

8.0 

1.779E+00 

1.102E+00 

8.134E-01 

6.763E-01 

5.796E-01 

5.118E-01 

1.023E-01 

9.0 

8.897E-01 

4.562E-01 

3.102E-01 

2.390E-01 

1.897E-01 

1.458E-01 

O.OOOE+OO 

10.0 

O.OOOE+00 

O.OOOE+OO 

O.OOOE+OO 

O.OOOE+OO 

O.OOOE+OO 

O.OOOE+OO 

O.OOOE+OO 

AE8MAX  MODEL 

B/BO 


L 

1.0 

2.0 

4.0 

6.0 

8.0 

10.0 

30.0 

1.0 

O.OOOE+OO 

O.OOOE+OO 

O.OOOE+OO 

O.OOOE+OO 

O.OOOE+OO 

O.OOOE+OO 

O.OOOE+OO 

2.0 

4.893E+00 

4.234E+00 

3.390E+00 

2.487E+00 

O.OOOE+OO 

O.OOOE+OO 

O.OOOE+OO 

3.0 

6.033E+00 

5.877E+00 

5.471E+00 

5.192E+00 

4.983E+00 

4.816E+00 

3.567E+00 

4.0 

6.505E+00 

6.405E+00 

6.221E+00 

6.086E+00 

5.979E+00 

5.913E+00 

5.470E+00 

5.0 

6.446E+00 

6.350E+00 

6.171E+00 

6.042E+00 

5.932E+00 

5.877E+00 

5.507E+00 

6.0 

5.999E+00 

5.903E+00 

5.724E+00 

5.599E+00 

5.489E+00 

5.438E+00 

5.101E+00 

7.0 

5.414E+00 

5.318E+00 

5.140E+00 

5.018E+00 

4.907E+00 

4.858E+00 

4.536E+00 

8.0 

4.462E+00 

4.366E+00 

4.187E+00 

4.066E+00 

3.955E+00 

3.907E+00 

3.593E+00 

9.0 

2.875E+00 

2.782E+00 

2.609E+00 

2.505E+00 

2.401E+00 

2.345E+00 

2.056E+00 

10.0 

6.982E-01 

6.288E-01 

4.899E-01 

4.125E-01 

3. 615E-01 

3.105E-01 

3.228E-02 
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4.  MODELS  OF  AURORAL  PRECEPITATION 


The  morphology  and  fluxes  of  precipitating  electrons  and  ions  in  the  Earth’s  auroral  region  is  an 
important  input  to  many  theoretical  calculations  and  is  critical  to  the  planning  of  auroral  imaging 
missions.  For  many  years,  the  DMSP  satellites  have  been  collecting  data  on  the  flux  of  energetic 
electrons  and  ions  over  the  poles.  These  data  have  been  accumulated  and  binned  as  a  function  of 
geomagnetic  coordinates  and  geophysical  parameters.  To  facilitate  the  use  of  these  binned  averages 
in  other  programs,  the  averages  have  been  fit  to  functional  forms.  These  models  result  in  compact 
and  efiScient  subroutines  by  which  fluxes  and  related  properties  can  be  evaluated. 

There  are  two  distinct  auroral  models  available.  Both  give  electron  and  ion  number  and  energy  flux 
as  a  function  of  corrected  geomagnetic  latitude  (A)  and  magnetic  local  time  (MLT).  They  differ  in 
that  one  set  of  models  assumes  that  the  variation  in  flux  at  a  specific  location  is  due  to  variation  in 
the  value  of  Kp.  The  other  set  assumes  that  the  variation  can  be  described  by  the  simultaneous 
specification  of  the  solar  wind  velocity  and  the  interplanetary  magnetic  field  component  B^.  These 
will  be  described  separately  although  their  use  and  output  is  basically  the  same. 

4.1  Models  Driven  by  Kp 

The  Kp  driven  models  consist  of  two  routines,  one  for  electrons  and  another  for  ions.  A  complete 
description  for  both  the  data  processing  and  modeling  for  the  electrons  can  be  found  in  Hardy,  et  at 
[1987].  The  data  processing  used  for  the  ion  model  generation  is  described  by  Hardy,  et  al.  [1989] 
and  the  model  generation  is  described  in  Hardy,  etal.  [1991].  However,  each  of  the  auroral  region 
“maps”  for  a  specific  value  of  Kp  is  treated  independently.  The  latitudinal  dependence  of  the  base  10 
logarithm  of  the  flux  at  each  discrete  MLT  is  examined  and  fit  to  an  Epstein  transition  function.  The 
first  order  Epstein  function,  which  was  used  for  the  electrons,  is  given  by 


e{X)  =  r  +  5j  (A  -  +  (Sj  -  Sj)  In^ 


1  -  S/5, 


This  complicated  expression  has  a  comparatively  simple  behavior,  since  e  is  equal  to  r  at  Ag  and  the 
curve  is  linear  with  slope  S,  and  S2  for  A  «  Aq  and  A  »Ao  respectively.  The  first  order  function  is 
used  for  electron  number  and  energy  flux  and  for  Hall  and  Pederson  conductivity.  A  second  order 
function  is  used  for  the  ion  number  and  energy  flux  because  the  curves  were  somewhat  broader  than 
the  electron  parameters  as  a  function  of  latitude.  The  second  order  Epstein  function  allows  two 
points  to  be  set,  namely  the  values  of  the  functions  at  Aq  and  Aj. 

The  coefficients  of  the  Epstein  functions  were  then  expanded  in  a  Fourier  series  for  each  map  to 
further  reduce  the  size  requirements  of  the  model  and  to  provide  another  degree  of  smoothing.  Each 
Kp  map  in  the  model  is  considered  separately  and  so  interpolation  in  this  dimension  is  not  available. 
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Two  FORTRAN  subroutines,  ‘elemod.f  and  Tonmod.f ,  are  available  to  evaluate  the  coefficients  of 
the  model  These  require  the  data  files  ELECOEF.DAT  and  IONCOEF.DAT.  Evaluation  of  a  full 
map  between  the  valid  latitudes  of  50°  and  90°  at  1°  increments  requires  about  8  seconds  with  ASCII 
output  and  about  4  seconds  without  output  on  a  Sparc  10  workstation.  Table  10  below  shows 
sample  output  for  ions  and  electron  fluxes  and  conductivities.  The  results  of  the  model  generation 
are  shown  in  Table  10  for  both  the  noon  and  midnight  slice  at  a  Kp  of  2. 


Table  10.  Sample  Output  from  Auroral  Kp  Models  for  Kp  of  2 

Midnight  Sector  Noon  Sector 


k 

Sp 

v 

Jn’ 

50.0 

0.00 

0.00 

l.OOOE+06 

l.OOOE+06 

52.0 

0.00 

0.00 

l.OOOE+06 

l.OOOE+06 

54.0 

0.00 

0.00 

l.OOOE+06 

l.OOOE+06 

56.0 

0.00 

0.00 

l.OOOE+06 

1.087E+06 

58.0 

0.00 

0.00 

2.112E+06 

2.982E+06 

60.0 

0.00 

0.00 

1.022E-K)7 

8.177E+06 

62.0 

4.77 

2.44 

4.812E+07 

2.233E-K)7 

64.0 

9.62 

5.40 

1.902E-K)8 

5.925E-K)7 

66.0 

9.93 

6.44 

3.767E+08 

1.320E+08 

68.0 

7.94 

5.78 

3.066E+08 

1.716E+08 

70.0 

5.49 

4.67 

1.795E+08 

1.360E+08 

72.0 

2.98 

3.48 

9.928E+07 

9.368E+07 

74.0 

0.55 

2.28 

5.446E+07 

6.308E+07 

76.0 

0.55 

1.07 

2.984E+07 

4.235E+07 

78.0 

0.55 

0.55 

1.635E+07 

2.842E+07 

80.0 

0.55 

0.55 

1.000E-K)7 

1.907E+07 

82.0 

0.55 

0.55 

l.OOOE+07 

1.280E+07 

84.0 

0.55 

0.55 

1.000E-K)7 

l.OOOE+07 

86.0 

0.55 

0.55 

l.OOOE+07 

l.OOOE+07 

88.0 

0.55 

0.55 

l.OOOE+07 

l.OOOE+07 

90.0 

0.55 

0.55 

l.OOOE+07 

l.OOOE+07 

Sh 


l.OOOE+04 

l.OOOE+04 

l.OOOE+04 

l.OOOE+04 

l.OOOE+04 

LOOOE+04 

2.248E+06 

5.017E+07 

5.981E+07 

3.412E+07 

1.630E+07 

7.186E-K)6 

3.100E-K)6 

L333E-K)6 

5.731E+05 

2.463E-K)5 

1.059E+05 

4.550E+04 

1.956E-K)4 

1.000E-K)4 

l.OOOE+04 


Je’ 


l.OOOE+04 

l.OOOE+04 

l.OOOE+04 

l.OOOE+04 

l.OOOE+04 

l.OOOE+04 

l.OOOE+04 

1.506E+04 

3.753E+04 

9.366E+04 

2.354E+05 

6.192E+05 

1.893E+06 

6.060E+06 

1.059E+07 

6.113E+06 

2.100E+06 

6.528E+05 

2.001E+05 

6.120E+04 

1.871E+04 


LOOOE+04 

l.OOOE+04 

l.OOOE+04 

LOOOE+04 

LOOOE+04 

1.163E+04 

2.898E+05 

2.450E+06 

4.602E+06 

3.166E+06 

1.636E+06 

8.074E+05 

3.959E+05 

1.939E+05 

9.498E+04 

4.652E+04 

2.279E+04 

1.116E+04 

l.OOOE+04 

l.OOOE+04 

l.OOOE+04 


0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

1.42 

0.41 

3.39 

1.14 

5.10 

1.84 

5.61 

2.32 

4.26 

2.14 

2.09 

1.34 

0.55 

0.55 

0.55 

0.55 

0.55 

0.55 

0.55 

0.55 

0.55 

0.55 

0.55 

0.55 

0.55 

0.55 

l.OOOE+06 

l.OOOE+06 

LOOOE+06 

l.OOOE+06 

l.OOOE+06 

LOOOE+06 

1.312E+06 

3.505E+06 

9.338E+06 

2.437E+07 

5.593E+07 

8.172E+07 

7.312E+07 

5.585E+07 

4.156E+07 

3.081E+07 

2.283E+07 

1.692E+07 

1.254E+07 

LOOOE+07 

LOOOE+07 


l.OOOE+06 

l.OOOE+06 

l.OOOE+06 

l.OOOE+06 

l.OOOE+06 

l.OOOE+06 

l.OOOE+06 

LOOOE+06 

l.OOOE+06 

3.020E+06 

9.583E+06 

3.028E+07 

9.264E+07 

2.326E+08 

3.022E+08 

2.018E+08 

1.099E+08 

5.794E+07 

3.040E+07 

1.594E+07 

l.OOOE+07 


l.OOOE+04 

l.OOOE+04 

l.OOOE+04 

LOOOE+04 

l.OOOE+04 

l.OOOE+04 

l.OOOE+04 

l.OOOE+04 

8.165E+04 

1.070E+06 

2,928E+06 

5.337E+06 

9.026E+06 

1.392E-K)7 

L459E+07 

8.070E+06 

3.316E+06 

1.287E+06 

4.956E+05 

1.906E+05 

7.330E+04 


4.2  Models  Driven  by  IMF  and  Solar  Wind 


Auroral  precipitation  models  in  the  other  set  are  driven  by  the  values  of  the  Interplanetary  Magnetic 
Field  (IMF)  and  the  Solar  Wind  Speed  (SWS).  Documentation  for  the  generation  of  these  models 
can  be  found  in  Brautigam,  etal  [1991].  They  differ  also  in  that  there  is  a  separate  set  of  data  files 
and  routines  for  the  conductivity.  The  fitting  was  done  with  the  same  functional  forms  except  that 
second  order  Epstein  functions  were  used  for  the  electron  fluxes  because  they  were  somewhat 
broader  than  in  the  Kp  versions.  Some  rough  comparisons  of  these  two  models  with  SSJ4  fluxes 
pass-by-pass  indicate  that  they  are  roughly  equivalent  in  accuracy,  however,  it  is  ultimately  the  user’ s 
responsibility  to  determine  which  set  of  models  best  serves  the  purpose.  This  is  to  say  that  neither 
the  Kp  or  the  IMF/S WS  model  is  recommended  over  the  other.  The  correspondence  between  map 
number  and  IMF/SWS  conditions  is  given  in  Table  1 1 . 


18 


Table  1 1 .  Map  Numbers  Corresponding  to  Various  IMF/SWS  Conditions 


B,\SWS 

346  km/s 

408  km/s 

485  km/s 

572  km/s 

677  km/s 

-4.5  nT 

map  1 

map  7 

map  13 

map  19 

map  25 

-2.2  nT 

map  2 

map  8 

map  14 

map  20 

map  26 

-0.7  nT 

map  3 

map  9 

map  15 

map  21 

map  27 

+0.7  nT 

map  4 

map  10 

map  16 

map  22 

map  28 

+2.2  nT 

map  5 

map  11 

map  17 

map  23 

map  29 

+4.5  nT 

map  6 

map  12 

map  18 

map  24 

map  30 

As  before,  the  models  are  in  the  form  of  subroutines  called  ionimf.f  eleimf.f,  and  conimf.f  and 
require  the  data  files  IONIMF.DAT,  ELEIMF.DAT,  and  CONIMF.DAT  respectively.  Time 
requirements  for  the  IMF/SWS  models  are  about  the  same  as  for  the  Kp  versions.  Table  12  below 
gives  the  expected  output  fi-om  map  15. 


Table  12.  Sample  output  from  Auroral  IMF/SWS  Models  for  Map  15 


Midnight  Sector 


Noon  Sector 


A 


Sh 


Sp 


Jn 


V 


Sp 


Je 


V  v 


50.0  0.50  0.25  3.732E+01  1.906E+02  l.OOOE+00  l.OOOE+00 
52.0  0.50  0.25  4,199E+02  1.447E-K)3  l.OOOE+00  l.OOOE+00 
54.0  0.50  0.25  4.725E'H)3  1.098E+04  8.895E-K)0  6.130E+00 
56.0  0.50  0.25  5.315E+04  8.332E-K)4  2.911E+02  1.314E+02 
58.0  0.50  0.25  5.968E+05  6.309E’K)5  9.434E+03  2.788E+03 
60.0  0.50  0.25  6.611E+06  4.708E-K)6  2.852E+05  5.464E+04 
62.0  3.54  2.97  6.660E+07  3.178E+07  5.794E+06  6.830E+05 
64.0  10.55  6.89  4.013E+08  1.324E+08  3.893E+07  2.574E-K)6 
66.0  16.01  9.36  7.314E+08  2.281E+08  7.301E+07  3.633E+06 
68.0  16.69  9.16  6.103E-H)8  2.509E+08  4.887E+07  3.003E+06 
70.012.99  7.58  4.250E+08  2.511E-H)8  2.204E+07  1.666E+06 
72.0  8.04  5.73  2.839E+08  2.387E+08  9.247E+06  8.109E+05 

74.0  2.90  3.83  1.708E+08  1.795E+08  3.841E+06  3.863E+05 

76.0  0.50  1.93  6.455E+07  6.499E+07  1.593E+06  1.834E+05 

78.0  0.50  0.25  1.124E+07  1.188E+07  6.606E-K)5  8.708E+04 

80.0  0.50  0.25  5.816E+06  7.644E+06  2.740E+05  5.438E+04 

82.0  0.50  0.25  5.816E+06  7.644E+06  2.326E+05  5.438E+04 

84.0  0.50  0.25  5.816E+06  7.644E+06  2.326E+05  5.438E+04 

86.0  0.50  0.25  5.816E+06  7.644E-K)6  2.326E+05  5.438E+04 

88.0  0.50  0.25  5.816E+06  7.644E+06  2.326E+05  5.438E+04 

90.0  0.50  0.25  5.816E+06  7.644E-K)6  2.326E+05  5.438E+04 


0.50 

0.25 

0.50 

0.25 

0.50 

0.25 

0.50 

0.25 

0.50 

0.25 

0,50 

0.25 

0.50 

0.25 

1.83 

0.25 

3.85 

0.72 

5.74 

1.20 

6.94 

1.67 

6.29 

2.09 

4.24 

2.23 

1.82 

1.83 

0.50 

1.13 

0.50 

0.38 

0.50 

0.25 

0.50 

0.25 

0.50 

0.25 

0.50 

0.25 

0.50 

0.25 

3.286E+01 

1.701E+02 

8.810E+02 

4.562E+03 

2.362E+04 

1.223E+05 

6.323E+05 

3.238E+06 

1.549E+07 

5.213E+07 

8.298E+07 

8.468E+07 

7.783E+07 

6.652E+07 

4.254E+07 

1.577E+07 

5.816E+06 

5.816E+06 

5.816E+06 

5.816E+06 

5.816E+06 


3.047E+01 

1.004E+02 

3.307E+02 

1.089E+03 

3.588E+03 

1.182E+04 

3.893E+04 

1.283E+05 

4.227E+05 

1.399E+06 

4.754E+06 

1.849E+07 

8.814E+07 

2.688E-K)8 

2.122E+08 

7.235E+07 

2.051E+07 

7.644E+06 

7.644E+06 

7.644E+06 

7.644E+06 


l.OOOE+00 

l.OOOE+00 

l.OOOE+00 

l.OOOE+00 

l.OOOE+00 

6.730E+00 

1.594E+02 

3.723E+03 

7.913E+04 

1.014E+06 

4.120E+06 

8.034E+06 

1.280E+07 

1.582E+07 

9.596E+06 

3.045E+06 

8.034E+05 

2.326E-K)5 

2.326E+05 

2.326E+05 

2.326E+05 


2.074E+02 

4.404E+02 

9.351E+02 

1.985E+03 

4.215E+03 

8.949E+03 

1.900E+04 

4.034E+04 

8.565E-K)4 

1.820E+05 

3.886E+05 

8.568E+05 

2.127E+06 

5.306E+06 

6.171E+06 

2.174E+06 

4.701E+05 

9.311E+04 

5.438E+04 

5.438E+04 

5.438E+04 
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5.  MAGNETIC  FIELD  MODELS;  CORRECTED  GEOMAGNETIC  COORDINATES 


The  Earth’s  magnetic  field  is  usually  modeled  as  the  sum  of  the  “main”  internal  magnetic  field,  which 
is,  for  the  most  part,  attributed  to  currents  within  the  earth’s  core,  and  an  “external”  field,  attributed 
to  ionospheric  and  magnetospheric  currents.  For  the  “main”  field  the  standard  reference  models  are 
the  IGRF/DGRF  models,  described  below.  The  Hilmer-Voight,  Tsyganenko-Usmanov,  and  Olson- 
Pfitzer  external  models  are  also  described  below. 

Corrected  Geomagnetic  (CGM)  Coordinates  are  based  upon  internal  magnetic  field  models,  and 
provide  a  convenient  tool  for  the  analysis  of  data  of  charged  particles  (electron,  protons  and  other 
ions)  moving  in  the  ionosphere  and  in  the  radiation  belts.  FORTRAN  routines  that  convert 
geographic  coordinates  to  CGM  and  the  inverse  are  also  described  below.  CGLALO  and 
CGLALO_INV  provide  table  interpolation  for  conversion  for  low  altitudes,  and 
SFC_CONVERT_GEO_COORD  provide  altitude  dependent  conversion  for  altitudes  0-7200  km. 
A  description  of  these  routines  is  provided  below. 


5.1  IGRF/DGRF  Internal  Magnetic  Field  Models 

The  International  Geomagnetic  Field  Reference  Field  models  represent  the  “main”  internal  magnetic 
field.  Since  the  internal  magnetic  field  changes  with  time,  the  model  is  updated  once  every  5  years 
by  the  International  Association  of  Geomagnetism  and  Aeronomy  (lAGA)  [Barton,  1996;  Langel, 
1991;  Langel,  etal,  1988].  Mathematically  the  model  is  represented  as  the  negative  gradient  of  a 
scalar  potential  V  represented  as  a  sum  of  spherical  harmonics  to  degree  and  order  10; 


V=aY. 

n=l 


10 


V  '•/ 


52  (^/’cos  m^)  P/’(sin  8) 

m-0 


(3) 


where  a  is  the  mean  earth  radius,  6371.2  km,  and  (r,  0,  (f))  are  the  spherical  coordinates  of  the 
position.  To  convert  fi'om  geodetic  to  spherical  coordinates,  lAGA  recommends  use  of  the 
International  Astronomical  Union  (lAU)  reference  ellipsoid,  with  equatorial  radius  6378.16  km  and 
flattening  ratio  1/298.25. 

The  models  have  been  derived  from  data  provided  by  observatories,  ships,  aircraft,  and  satellites. 
Daily  and  annual  averages  of  the  measurements  are  used  in  least  squares  fits  to  compute  the  spherical 
harmonic  coefficients.  If  possible,  daily  values  used  for  computing  the  fits  are  limited  to  magnetically 
quiet  days,  to  avoid  perturbations  by  external  currents.  This  problem  is  particularly  troublesome  for 
satellite  data.  Many  satellites  pass  through  regions  of  external  currents,  for  which  the  spherical 
harmonic  representation  is  not  valid  (curl  of  B  is  not  zero).  Principally  there  exist  field-aligned 
currents,  running  parallel  to  field  lines  between  the  magnetosphere  and  the  high-latitude  ionosphere. 


The  field  aligned  currents  perturb  mainly  the  direction,  but  not  the  magnitude  of  the  field.  Therefore 
only  the  magnitude  is  used  in  the  analysis.  Langel  and  Estes  [1985]  found  that  eastward  flowing 
auroral  electrojets  underneath  satellite  MAGS  AT  produced  asymmetries  inlocaltime,  with  significant 
differences  between  dawn  and  dusk  measurements.  Data  for  the  subsequent  analysis  were  therefore 
carefully  selected,  giving  precedence  to  dawn,  where  the  effect  was  minimal.  Where  necessary,  data 
fi'om  dusk  were  used  with  latitude-dependent  corrections. 

The  existing  models  consists  of  a  series  of  DGRF  (Definitive  Geomagnetic  Reference  Field)  models 
at  5-year  intervals  from  1945  through  1990,  and  an  IGRF  model  for  1995,  where  the  latter  includes 
coefficients  for  the  predicted  secular  variation  for  the  years  1995  thru  2000.  Each  DGRF  model  is 
specified  by  a  set  of  coefficients  h„”.  Values  of  the  coefficients  between  successive  DGRF  epochs 
are  obtained  by  linear  interpolation.  The  IGRF  models  include  the  latter  and  estimates  of  their  first 
time  derviatives.  The  DGRF  models  are  expected  not  to  change,  since  no  new  data  is  likely  to  be 
available  for  them.  For  the  period  between  1990  and  1995  interpolation  is  performed  between 
DGRF90  and  the  set  of  coefficients  hn”  from  IGRF95. 

In  the  year  2000,  the  IGRF95  model  will  be  replaced  by  a  DGRF95  model  based  on  new 
observational  data,  and  a  new  predictive  model  IGRF  2000  will  become  available,  with  an  estimated 
secular  variation  for  the  following  5  years. 

To  take  into  account  the  effects  of  the  external  fields,  it  is  customary  to  augment  the  DGRF  and 
IGRF  models  with  external  field  models,  some  of  which  are  described  below.  The  IGRF  models  also 
form  the  basis  of  computations  of  Corrected  Geomagnetic  Coordinate  routines  such  as  CGLALO  and 
SFC_CONVERT_GEO_COORD. 

Several  codes  are  available  to  compute  both  the  magnetic  field  and  derived  properties,  such  as  L- 
shell.  We  use  here  the  NASA  package,  which  employs  the  inverse  coordinate  method  of  Kluge 
[1972].  Table  13  presents  a  sample  output  for  1  Jan,  0  hr  UT,  in  the  indicated  years.  The  total 
magnetic  field,  and  its  northward,  eastward  and  downward  components  are  given  in  Gauss.  The  L- 
computation  routine  returns  a  code  (ICODE)  which  is  1  to  indicate  successful  computation  by  the 
normal  method  of  numerical  integration,  2  to  indicate  failure  because  the  conjugate  point  is  below 
the  earth’s  surface  (the  value  output  is  meaningless),  and  3  to  indicate  that  L  was  computed  by  an 
approximate  method,  applicable  at  very  large  L. 

On  a  40  MHZ  Sparc  IPX  workstation,  the  average  times,  per  10000  calls,  for  the  various 
computations  are; 

Model  coefficients  102  sec 

Magnetic  field  vector  0.7  sec 

L  7.5  sec 
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Table  13.  Sample  IGRF90  Model  Output 


YEAR 

ALTITUDE 

LAT 

LON 

B (NORTH) 

B  (EAST) 

B  (DOWN) 

B (TOTAL) 

L 

I  CODE 

1946 

800 

40 

0 

0.168976 

-0.027447 

0.254116 

0.306400 

1.7909 

1 

1956 

800 

40 

0 

0.171061 

-0.024426 

0.254050 

0.307246 

1.7773 

1 

1966 

800 

40 

0 

0.173115 

-0.021878 

0.253273 

0,307563 

1.7608 

1 

1976 

800 

40 

0 

0.175535 

-0.019309 

0,252736 

0.308320 

1.7381 

1 

1986 

800 

40 

0 

0.177165 

-0.015420 

0.251816 

0.308280 

1.7198 

2 

1996 

800 

40 

0 

0.178122 

-0.012116 

0.251588 

0.308498 

1.7068 

2 

2006 

0 

40 

0 

0.253320 

-0.010634 

0.363182 

0.442928 

1.4601 

2 

2006 

5000 

40 

0 

0.042920 

-0.005833 

0.063948 

0.077237 

2.9277 

1 

2006 

10000 

40 

0 

0.013862 

-0.002348 

0.022052 

0.026153 

4.3749 

1 

2006 

15000 

40 

0 

0.006101 

-0.001142 

0.010101 

0.011856 

5.8417 

1 

2006 

20000 

40 

0 

0.003205 

-0.000635 

0.005443 

0.006348 

7.3239 

1 

2006 

25000 

40 

0 

0.001887 

-0.000388 

0.003262 

0.003788 

8.7915 

1 

2006 

800 

15 

0 

0.221757 

-0.019374 

0.046358 

0.227378 

1.1294 

1 

2006 

800 

25 

0 

0.216404 

-0.015066 

0.137736 

0.256961 

1.2289 

1 

2006 

800 

35 

0 

0.193437 

-0.012634 

0.217962 

0.291694 

1.4844 

1 

2006 

800 

45 

0 

0.161603 

-0.011998 

0.280723 

0.324137 

2.0204 

2 

2006 

800 

55 

0 

0.127604 

-0.012523 

0.326832 

0.351082 

3.0926 

2 

2006 

800 

65 

0 

0.094509 

-0.013060 

0.360054 

0.372480 

5.6217 

1 

2006 

800 

40 

50 

0.182531 

0.010100 

0.281426 

0.335590 

1.6364 

1 

2006 

800 

40 

100 

0.196909 

-0.003280 

0.311539 

0.368566 

1.5330 

1 

2006 

800 

40 

150 

0.193310 

-0.010558 

0.255774 

0.320781 

1.5213 

1 

2006 

800 

40 

200 

0.166680 

0.037719 

0.260741 

0.311754 

1.7902 

1 

2006 

800 

40 

250 

0.145364 

0.029618 

0.337651 

0.368804 

2.4993 

1 

2006 

800 

40 

300 

0.142618 

-0.042709 

0.318770 

0.351822 

2.6275 

1 

5.2  Hilmer-Voigt  Magnetospheric  Magnetic  Field  Model 


Unlike  the  internal  field,  there  is  no  standard  reference  model  for  the  external  field.  Several 
models  have  been  developed  with  varying  degrees  of  sophistication.  Among  the  most  elegant  is 
the  Hilmer-Voigt  model  [Hilmer  and  Voigt,  1995;  Hilmer,  1989].  This  model  seeks  to  provide 
realistic,  yet  mathematically  tractable,  representations  of  known  essential  features  of  the 
magnetosphere.  The  condition  V  •  B  “  0  is  satisfied  throughout.  Contributions  due  to 
magnetopause  currents,  ring  currents,  and  cross-tail  currents  are  included.  The  magnetopause  field 
in  the  Hilmer-Voigt  model  "shields"  the  dominant  dipole  portion  of  the  internal  field  from  the 
solar  wind  by  canceling  the  dipole  field’s  normal  component  along  the  magnetopause  [Voigt, 
1981].  The  magnetopause  consist  of  a  hemisphere  facing  the  sun  with  a  cylinder  of  equal  radius 
extending  anti-sunward.  The  ring  currents  consist  of  eastward  and  westward  flowing  currents, 
cylindrically  symmetric  with  respect  to  the  earth.  They  have  a  simple  mathematical  form,  first 
used  by  Tsyganenko  and  Usmanov  [1982],  whose  model,  however,  has  only  a  westward  flowing 
current.  The  eastward  current,  at  distances  less  than  3-4  Rg,  is  inferred  from  observations  of  a 
maximum  in  the  plasma  pressure,  which  in  equilibrium  satisfies  the  condition  J  x  B  =  VP.  The 
maximum  in  P  forces  VP  to  change  sign,  which  means  that  the  current  density  must  change  sign. 
The  cross-tail  currents,  embedded  in  the  night-side  plasma  sheet,  are  represented  by  current 
filaments,  flowing  infinitely  from  dawn  to  dusk,  in  a  series  of  16  plane  segments.  The  segments 
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are  positioned  and  oriented  to  simulate  a  sheet  with  a  specified  near-earth  edge  which  is  in  the 
earth’s  magnetic  equator.  When  the  earth’s  dipole  is  not  tilted  with  respect  to  the  earth-sun  line, 
the  entire  sheet  remains  on  the  equator.  When  the  dipole  is  tilted,  the  sheet  curves  toward  an 
orientation  which,  at  large  distances  from  the  earth,  is  parallel  to  the  GSM  equator  (the  plane 
through  the  earth-sun  line  that  is  normal  to  the  plane  containing  the  sun  and  the  dipole). 

The  model  is  driven  by  four  parameters:  dipole  tilt,  magnetospheric  standoff  distance  (noon 
equatorial  intersection  of  the  magnetopause),  geomagnetic  activity  index  Djt,  and  the  noon 
equatorward  boundary  of  the  diffuse  aurora.  The  standoff  distance  and  dipole  tilt  determine  the 
magnetopause  field.  The  ring  and  cross-tail  currents  are  described  by  structural  and  strength 
parameters  which  are  determined  by  a  procedure  that  forces  the  model  field  to  reproduce  and 
known  AB  patterns  in  the  noon-midnight  meridian,  while  also  satisfying  the  requirement  that  the 
cross-tail  inner  edge  maps  to  the  equatorward  boundary  of  diffuse  aurora.  In  practice  the  strength 
and  structural  parameters  are  prestored  on  file  for  specific  configurations.  At  run  time  the 
configuration  closest  to  the  specified  set  of  inputs  (tilt,  standoff,  Djt,  and  equatorward  boundary) 
is  chosen. 

An  option  exists  to  simulate  substorm  diversion  of  the  cross-tail  currents  near  the  noon-midnight 
meridian,  by  adding  an  eastward  current  to  the  normal  westward  current.  One  strength  parameter 
and  two  structure  parameters  are  determined  to  reproduce  a  specific  midnight  geosynchronous 
recovery  of  and  to  maintain  positive  B^  over  the  entire  equator. 

The  sample  output  shown  in  Table  14  was  generated  at  8  locations,  for  two  configurations  defined 
by  the  parameters  shown:  standoff  distance,  tilt,  Djj,  and  equatorward  edge  (EQEDGE).  The  8 
locations  are  specified  in  the  GSM  coordinate  system.  On  a  40MHz  Sparc  workstation,  the 
Hilmer-Voigt  model  takes  0.3  seconds  for  initialization,  which  must  take  place  whenever  any  of 
the  model  input  parameters  changes.  The  inititialization  includes  setting  up  the  structural  and 
strength  parameters,  based  on  the  set  of  input  parameters,  and  performance  of  numerical  tasks 
which  depend  on  these  parameters.  Initialization  -f  10000  calls  to  compute  point  model  magnetic 
field  vectors  requires  3.2  sec. 
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Table  14.  Sample  Hilmer-Voigt  Output 


Standoff  = 

10.00  RE 

Tilt  =  10.00 

Deg  DST 

= 

-20.0 

nT  EQEDGE  =  68.0 

X  (RE) 

Y  (RE) 

Z  (RE) 

BX  (nT) 

BY 

(nT) 

BZ  (nT) 

BTOT  (nT) 

2. 

0. 

0. 

-1334.7 

0.0 

3765.0 

3994.6 

0. 

2. 

0. 

664.9 

-0.3 

3761.8 

3820.1 

0. 

0. 

2. 

668.2 

0.0 

-7584.7 

7614.1 

6. 

0. 

0. 

-47.9 

0.0 

143.5 

151.3 

0. 

6. 

0. 

24.5 

-0.8 

132.5 

134.7 

0. 

0. 

6. 

33.6 

0.0 

-302.6 

304.5 

-2. 

0. 

0. 

-1336.4 

0.0 

3758.4 

3988.9 

-6. 

0. 

0. 

-54.0 

0.0 

121.3 

132.8 

Standoff  = 

7.00  RE 

Tilt  =  10.00 

Deg  DST 

=  ” 

•200.0 

nT  EQEDGE  =  55.0 

X  (RE) 

Y  (RE) 

Z  (RE) 

BX  (nT) 

BY 

(nT) 

BZ  (nT) 

BTOT  (nT) 

2. 

0. 

0. 

-1362.9 

0.0 

3661.8 

3907.2 

0. 

2. 

0. 

645.7 

-2.1 

3638.7 

3695.5 

0. 

0. 

2. 

663.8 

0.0 

-7796.6 

7824.8 

6 . 

0. 

0. 

-11.2 

0.0 

301.6 

301.8 

0. 

6. 

0. 

30.8 

-4.6 

165.8 

168.7 

0. 

0. 

6. 

71.1 

0.0 

-332.9 

340.4 

-2. 

0. 

0. 

-1374.7 

0.0 

3622.4 

3874.5 

-6. 

0. 

0. 

-71.8 

0.0 

144.8 

161.6 

Deg 


5.3  Olson-Pfitzer  Quiet  Tilt-Dependent  External  Field  Model 

The  Olson-Pfitzer  tilt  dependent  external  field  model  [Olson  and  Ffitzer,  1977]  includes  the  field 
contributions  from  the  magnetopause,  tail,  and  ring  currents.  These  contributions  are  formed  and 
oriented  by  the  solar  wind,  and  are  thus  tilted  as  a  function  of  season  and  time,  independent  of 
Earth-dipole  axis  coordinates.  Tilt  is  defined  as  the  complement  of  the  Sun-dipole  angle.  The  total 
external  field  is  obtained  by  vectorially  adding  this  component  to  the  internal  field,  which  is 
defined  by  a  standard  model  such  as  IGRF  1985  [Barraclough,  1987].  The  Olson-Pfitzer  model 
revises  an  earlier  tilt-averaged  model  that  did  not  adequately  represent  the  zero  tilt  case, 
particularly  at  large  distances  on  the  midnight  equator. 


This  model  is  valid  for  all  tilts  of  the  earth’s  dipole.  There  is  no  dependence  on  magnetic  activity 
-  the  model  is  an  average  for  quiet  conditions  (Kp  =  0,  1).  The  code  has  been  optimized  for  the 
near  earth  region,  and  is  valid  to  15  Rg.  Beyond  that,  the  field  diverges  rapidly,  and  a  template 
sets  the  field  to  zero.  To  improve  computational  speed,  the  external  field  is  set  to  zero  below  2 
Re,  since  the  external  field  contribution  to  the  total  field  is  negligible  in  this  region. 


Magnetopause  currents  and  their  resulting  contributions  to  the  field  are  computed  numerically; 
ring  and  cross-tail  currents  are  constructed  from  a  conceptual  system  of  wires.  Parameters 
specifying  the  wires  (position,  orientation,  size,  and  shape),  their  currents,  and  the  magnetopause 
shape  are  adjusted  so  that  the  resulting  total  external  field  fits  observations  from  the  inner 
magnetosphere  and  the  tail.  The  resulting  field  is  then  fit  by  linear  least  squares  to  an  analytic 
series  [Bass  and  Jordan,  1990]. 
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The  sample  runs  summarized  in  Table  15  were  generated  using  code  incorporating  the  IGRF 1985 
internal  field  model.  These  examples  show  the  model’s  seasonal  variation,  as  well  as  the  variation 
with  hour  (UT),  radius  (in  km),  and  latitude.  Program  output  consists  of  the  components  of  the 
model  magnetic  field  in  geocentric  coordinates,  along  with  the  total  field  (in  nanoTesla) .  The  time 
required  for  the  external  field  calculation  is  small:  10,000  total  field  calculations  require  57 
seconds  on  a  40  MHz  Sparc  IPX  workstation,  while  51  seconds  are  needed  to  compute  only  the 
internal  field. 


Table  15.  Olson-Pfitzer  Tilt  Dependent  External  Field  Model  Sample  Runs 


Year 

Day 

UT 

Radius 

Lat 

Lon 

Bx 

By 

Btot 

1994 

1 

0 

50000 

40 

0 

-112.153 

-8.327 

-40.400 

119.50 

1994 

50 

0 

50000 

40 

0 

-116.544 

-9.165 

-38.527 

123.09 

1994 

100 

0 

50000 

40 

0 

-122.376 

-8.504 

-35.006 

127.57 

1994 

150 

0 

50000 

40 

0 

-126.346 

-8.007 

-31.660 

130.50 

1994 

200 

0 

50000 

40 

0 

-126.256 

-8.738 

-31.913 

130.52 

1994 

250 

0 

50000 

40 

0 

-121.850 

-8.262 

-35.283 

127.12 

1994 

300 

0 

50000 

40 

0 

-115.761 

-7.586 

-38.585 

122.26 

1994 

100 

0 

50000 

40 

0 

-122.376 

-8.504 

-35.006 

127.57 

1994 

100 

4 

50000 

40 

0 

-111.975 

4.698 

-24.178 

114.65 

1994 

100 

8 

50000 

40 

0 

-95.599 

0.748 

-7.170 

95.87 

1994 

100 

12 

50000 

40 

0 

-90.374 

-15.514 

3.603 

91.77 

1994 

100 

16 

50000 

40 

0 

-97.813 

-30.260 

-9.118 

102.79 

1994 

100 

20 

50000 

40 

0 

-116.311 

-26.549 

-31.476 

123.39 

1994 

100 

24 

50000 

40 

0 

-122.480 

-8.482 

-34.927 

127.64 

1994 

100 

0 

20000 

40 

0 

-1429.696 

-128.736 

-256.545 

1458.22 

1994 

100 

0 

30000 

40 

0 

-444.255 

-38.485 

-106.030 

458.35 

1994 

100 

0 

40000 

40 

0 

-207.373 

-16.408 

-56.742 

215.62 

1994 

100 

0 

50000 

40 

0 

-122.376 

-8.504 

-35.006 

127.57 

1994 

100 

0 

60000 

40 

0 

-84 . 982 

-4.962 

-23.412 

88.29 

1994 

100 

0 

70000 

40 

0 

-66 .576 

-3.885 

-14.995 

68.35 

1994 

100 

0 

80000 

40 

0 

-55.129 

-4.767 

-6.257 

55.69 

1994 

100 

0 

50000 

-80 

0 

62 .894 

-12.968 

-122.470 

138.28 

1994 

100 

0 

50000 

-60 

0 

106.151 

-10.897 

-89.105 

139.02 

1994 

100 

0 

50000 

-40 

0 

108.226 

-9.605 

-30.936 

112.97 

1994 

100 

0 

50000 

-20 

0 

64.434 

-8.843 

25.109 

69.72 

1994 

100 

0 

50000 

0 

0 

-10.325 

-8.356 

47.366 

49.19 

1994 

100 

0 

50000 

20 

0 

-83.513 

-8.179 

21.993 

86.75 

1994 

100 

0 

50000 

40 

0 

-122.376 

-8.504 

-35.006 

127.57 

1994 

100 

0 

50000 

60 

0 

-112.236 

-9.804 

-91.097 

144.88 

1994 

100 

0 

50000 

80 

0 

-59.826 

-12.180 

-119.158 

133.89 
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5.4  Tsyganenko-Usmanov  External  Field  Model 

The  Tsyganenko-Usmanov  external  field  model  [Tsyganenko  and  Usmanov,  1982]  is  expressed  as 
the  sum  of  contributions  firom  ring  currents,  cross-tail  currents,  and  magnetopause  currents.  No 
attempt  was  made  to  provide  an  accurate  model  of  the  near-Earth  region,  due  to  limited  data 
coverage,  but  this  model  provides  reasonable  results  fi'om  5  to  20  R^.  This  model  is  the  forerunner 
of  the  Tsyganenko  model  [Tsyganenko,  1987], 

The  ring  and  tail  contributions  are  expressed  analytically  by  functions  that  are  curls  of  vector 
potentials;  therefore,  their  divergences  vanish  identically.  The  model  employs  a  mathematically  simple 
azimuthally  symmetric  ring  current  circulating  around  the  dipole.  The  ring  current  flow  direction  is 
westward,  independent  of  distance  from  the  earth.  The  tail  current  is  represented  as  flowing  mainly 
from  dawn  to  dusk  on  a  flat  “sheet”.  The  contribution  of  the  magnetopause  current  is  represented 
by  tilt-dependent  products  of  polynomials  and  decaying  exponentials  and  Jordan,  1990]. 

The  database  for  the  Tsyganenko-Usmanov  model  consists  of  data  from  the  Mead-Fairfield  model, 
and  data  from  HEOS-1  and  HEOS-2.  The  Mead-Fairfield  data  base  consists  of  averaged 
measurements  fi'om  magnetometers  on  board  Explorer  33,  Explorer  34,  Explorer  41,  and  Explorer 
43.  These  data  were  averaged  over  half-Earth  radii  -  typically  10  to  15  minutes.  The  Tsyganenko- 
Usmanov  data  set  has  a  gap  in  the  near-Earth  region,  making  it  unreliable  inside  of  5  Rg.  There  is 
also  minimal  data  coverage  beyond  20  R^,  defining  the  model’s  practical  limit. 


5.5  Corrected  Geomagnetic  Coordinate  Codes:  CGLALO,  GGLALO_INVERSE, 
SFC_CONVERT_GEO_COORD 

Ionospheric  processes  and  the  motion  of  charged  particles,  both  in  the  ionosphere  and  the  radiation 
belts,  are  determined  by  Earth’s  magnetic  field.  In  the  early  period  of  space  exploration,  the  use  of 
Geomagnetic  Coordinates,  a  spherical  coordinate  system  aligned  with  an  earth  centered  dipole,  was 
suitable  for  the  analysis  of  rocket  and  satellite  based  observations  of  the  motion  of  charged  particles 
and  ionospheric  phenomena.  As  improvements  in  observational  techniques  were  developed,  and  more 
precise  measurements  became  possible,  it  became  necessary  to  develop  coordinates  that  more  closely 
reflected  the  actual  magnetic  field  of  the  Earth.  These  Corrected  Geomagnetic  Coordinates  (CGM) 
systems  were  based  upon  magnetic  field  line  tracing  using  magnetic  field  models,  in  which  the 
magnetic  field  is  represented  as  the  negative  gradient  of  a  magnetic  scalar  potential. 

The  original  procedure  to  compute  CGM  coordinates  was  as  follows:  A  real  field  line  from  the 
Earth’s  surface  is  traced  to  the  equator  of  the  centered  dipole.  This  point  is  defined  to  be  the 
equivalent  to  a  line  trace  along  a  centered  dipole  field.  The  latitude  and  longitude  of  the  landing  point 
in  centered  dipole  coordinates  are  the  desired  “corrected”  geomagnetic  coordinates. 
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CGM  coordinates,  unlike  geographic  coordinates  or  centered  dipole  magnetic  coordinates,  are 
distorted,  the  distortion  being  a  reflection  of  the  deviation  of  the  Earth’s  magnetic  field  from  a 
centered  dipole  field.  An  Earth  globe  based  upon  CGM  coordinates  will  distort  distances  and  angles 
between  pairs  of  points  compared  to  one  based  upon  geographic  coordinates.  This  distortion  will  be 
particularly  evident  in  the  South  Atlantic  Anomaly  (SAA),  where  the  actual  (offset)  location  of  the 
dipole  component  of  Earth’s  magnetic  field  is  closest  to  Earth’s  surface. 

The  CGM  coordinates  for  a  given  magnetic  field  model  were  provided  to  users  in  the  form  of  a  table 
for  a  uniform  geographic  coordinate  grid  (uniform  spacing  in  latitude  and  longitude).  Some  of  these 
tables  were  incomplete,  and  excluded  low  (CGM)  latitudes.  CGLALO  incorporates  such  a  table, 
using  four  point  interpolation  to  compute  the  CGM  coordinates  except  for  the  vicinity  of  the  CGM 
poles  where  a  three  point  interpolation  method  is  used,  with  the  CGM  pole  being  one  of  the  three 
points.  The  CGLALO_INVERSE  routine  computes  geographic  coordinates  from  CGM  coordinates 
by  using  a  search  routine,  and  then  performing  a  Newton-Ralphson  procedure  using  the  CGLALO 
routine.  The  current  version  of  these  two  routines  CGLAL095  and  CGL95INV  are  based  upon  the 
IGRF  1995  Magnetic  Field  Model.  To  alleviate  (and  smooth)  distortions  in  the  vicinity  of  the 
magnetic  equator,  the  IGRF  1995  tables  used  to  generate  CGLALO  were  modified  as  follows: 

The  magnetic  dip  equator  is  the  set  of  points  on  the  Earth  reference  sphere  at  which  the  actual 
magnetic  field  lines  (as  defined  by  the  magnetic  field  model)  are  tangent  to  the  sphere.  Note  that  the 
dip  equator  does  not  actually  lie  in  a  plane,  but  is  taken  to  be  the  reference  plane  (equator)  for  CGM 
coordinates.  From  the  original  tables,  the  entries  within  a  band  of  ~15  degrees  around  the  magnetic 
dip  equator  were  deleted,  and  replaced  by  values  obtained  using  a  spline  fit,  with  the  added  constraint 
that  the  spline  curves  pass  through  the  coordinates  of  the  dip  equator. 

The  SFC_CONVERT_GEO_COORD  routine  computes  Altitude  Adjusted  Corrected  Geomagnetic 
Coordinates  fi-om  Geographic  Coordinates  and  altitude,  and,  where  they  exist,  the  inverse  from  CGM 
coordinates  and  altitude.  In  this  enhancement  of  the  CGM  coordinate  system  concept,  the  CGM 
coordinates  of  an  actual  field  line  are  constant  for  the  portion  of  the  field  line  north  or  south  of  the 
dip  equator  (at  0  km  altitude).  The  CGM  coordinates  are  the  landing  points  of  the  field  line,  as 
defined  by  the  CGLALO  routine. 

The  original  version  was  developed  by  Baker  and  Wing  [1989]  where  the  corrected  geomagnetic 
coordinates  (and  the  corresponding  inverses)  are  computed  by  evaluation  of  functions  for  the  X,  Y, 
and  Z  components  of  a  unit  vector  obtained  from  a  spherical  harmonic  expansion.  Since  they  were 
primarily  interested  with  representing  higher  CGM  latitudes,  the  South  Atlantic  Anomaly  and  the 
CGM  equatorial  region  were  not  well  represented.  Subsequent  versions  are  based  upon  an  improved 
method  described  in  Bhavnani  and  Hein  [1994].  These  versions  were  limited  to  an  altitude  range  of 
0-2000  km.  The  current  version  is  based  upon  additional  improvements  described  in  Hein  and 
Bhavnani  [1996]  has  an  altitude  range  of  0-7200  km.  Versions  of  the  current  routine,  and  the 
CGLALO  routines  based  upon  the  DGRF  1975, 1980, 1985, 1990andIGRF  1995  models  have  been 
prepared  for  the  research  community. 
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Since  the  altitude  adjusted  CGM  coordinates  along  a  field  line  trace  are  (by  definition)  constant  on 
a  CGM  hemisphere,  and  differ  only  by  the  sign  of  CGM  latitude  as  they  cross  hemispheres,  it  is 
possible  to  use  the  current  version  (and,  with  less  accuracy,  previous  versions)  for  field  line  tracing. 
A  description  of  such  usage  and  comparison  with  precision  field  line  tracing  routines  are  provided  in 
Hein  cmdBhamani  [1996],  Table  16  provides  a  comparison  of  this  usage  with  calculations  using  a 
precision  magnetic  field  line  tracing  routine. 


Table  16. 

Comparison  of  Field  Line  Trace  application  of  SFC_CONVERT_GEO_COORDINATE  Routine 

with  Precise  Magnetic  Field  Line  Trace. 


start  Location  for 

Trace 

SFC 

CONVERT.  ' 

output 

PRECISE 

F.L.  Trace 

output 

Altitude 

G.Lat. 

G.  Loner 

Alt 

G.Lat. 

G .  Loner  - 

G.  Lat. 

G.  Lon. 

Diff. 

[km] 

[deg] 

[deg] 

[km] 

[deg] 

[deg] 

[deg] 

[deg] 

[deg] 

7200.0 

30.000 

.000 

.0 

55.938 

355.040 

55.288 

355.300 

.667 

7200.0 

30.000 

30.000 

.0 

55.906 

29.272 

55.750 

29.253 

.157 

7200.0 

30.000 

60.000 

.0 

55.358 

60.747 

55.298 

60.693 

.067 

7200.0 

30.000 

90.000 

.0 

55.033 

89.938 

55.098 

89.987 

.071 

7200.0 

30.000 

120.000 

.0 

55.769 

118.530 

55.765 

118.512 

.010 

7200.0 

30.000 

150.000 

.0 

57.012 

150.055 

57.089 

150.150 

.093 

7200.0 

30.000 

180.000 

.0 

56.040 

185.768 

56.099 

185.812 

.064 

7200.0 

30.000 

210.000 

.0 

52.549 

217.958 

52.534 

217.960 

.014 

7200.0 

30.000 

230.000 

.0 

50.208 

236.728 

50.184 

236.657 

.052 

7200.0 

30.000 

270.000 

.0 

47.526 

270.203 

47.519 

270.227 

.018 

7200.0 

30.000 

300.000 

.0 

48.448 

294.408 

48.655 

294.484 

.213 

7200.0 

30.000 

330.000 

.0 

52.667 

321.918 

52.334 

322.154 

.363 

7200.0 

60.000 

.000 

.0 

70.789 

354.165 

70.685 

354.462 

.143 

7200.0 

60.000 

30.000 

.0 

71.284 

28.946 

71.073 

29.236 

.230 

7200.0 

60.000 

60.000 

.0 

70.585 

61.687 

70.487 

61.670 

.098 

7200.0 

60.000 

90.000 

.0 

69.963 

90.753 

69.980 

90.893 

.051 

7200.0 

60.000 

120.000 

.0 

70.386 

119.247 

70.295 

119.252 

.091 

7200.0 

60.000 

150.000 

.0 

71.149 

150.897 

71.110 

150.630 

.095 

7200.0 

60.000 

180.000 

.0 

70.755 

185.718 

70.817 

185.635 

.067 

7200.0 

60.000 

210.000 

.0 

68.741 

217.917 

68.851 

218.055 

.121 

7200.0 

60.000 

230.000 

.0 

67.326 

236.741 

67.383 

236.624 

.073 

7200.0 

60.000 

270.000 

.0 

65.770 

269.499 

65.949 

269.466 

.179 

7200.0 

60.000 

300.000 

.0 

66.775 

293.636 

66.825 

293.862 

.102 

7200.0 

60.000 

330.000 

.0 

68.985 

321.455 

68.907 

321.661 

.107 

.0 

50.000 

.000 

3600.0 

33.782 

1.955 

33.817 

1.988 

.045 

.0 

50.000 

30.000 

3600.0 

34.239 

29.708 

34.166 

29.716 

.073 

.0 

50.000 

60.000 

3600.0 

35.133 

59.033 

35.172 

59.055 

.043 

.0 

50.000 

90.000 

3600.0 

35.657 

89.797 

35.587 

89.793 

.071 

.0 

50.000 

120.000 

3600.0 

34.549 

121.489 

34.513 

121:449 

.049 

.0 

50.000 

150.000 

3600.0 

31.833 

150,677 

31.792 

150.658 

.044 

.0 

50.000 

180.000 

3600.0 

32.461 

177.091 

32.609 

177.099 

.148 

.0 

50.000 

210.000 

3600.0 

36.410 

205.000 

36.533 

205.135 

.164 

.0 

50.000 

230.000 

3600.0 

39.188 

225.403 

39.205 

225.518 

.091 

.0 

50.000 

270.000 

3600.0 

42.382 

270.049 

42.267 

270.066 

.116 

.0 

50.000 

300.000 

3600.0 

40.750 

304.276 

40.645 

304.160 

.137 

.0 

50.000 

330.000 

3600.0 

36.771 

334.628 

36.483 

334.593 

.289 

The  Diff.  Field  represents  the  difference  along  the  arc  of  a  great  circle  passing 
through  the  two  points. 
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6.  ATMOSPHERE/THERMOSPHERE  MODELS 


The  models  described  in  this  section  provide  estimates  of  temperature,  density  and  composition 
(neutral  atoms  and  molecules)  as  a  function  of  local  time,  day  of  year,  altitude,  geographic  location, 
and  geophysical  conditions  (solar  and  geomagnetic  activity).  Although  the  ionosphere  lies  within  the 
model  ranges,  the  free  electron  and  ion  content  are  not  part  of  these  models,  but  are  described  in  the 
section  on  ionospheric  models.  References  and  sample  test  output  for  these  models  is  provided. 


6.1  MSIS  Thermosphere  Models 

These  models  are  empirical  models  based  on  mass  spectrometer  and  incoherent  scatter  data  -  hence, 
the  acronym  MSIS.  Mass  spectrometer  data,  taken  from  instruments  on  board  rockets  and  satellites, 
provides  measurements  of  the  number  densities  of  the  major  constituents  in  the  altitude  range  100- 
700  km.  Temperature  data  were  taken  from  these  satellites,  as  well  as  from  rocket  measurements 
between  85  and  300  km,  and  from  ground  based  measurements  of  incoherent  scatter  at  low  altitudes 
(110-400  km)  and  at  exospheric  (above  500  km)  altitudes. 

The  MSIS  thermospheric  models  generate  temperature  and  densities  to  an  altitude  of  1 000  km,  using 
analytic  vertical  temperature  profiles.  These  temperature  profiles  depend  on  geographic  coordinates, 
local  time,  day  of  year,  and  geophysical  conditions  (solar  flux  and  geomagnetic  activity).  Vertical 
profiles  of  the  densities  are  obtained  from  the  temperature  profiles  as  solutions  of  differential 
equations  for  diffusive  and  mixing  equilibrium.  In  the  mixing  profiles,  the  densities  of  all  constituents 
decay  at  the  same  rate  with  altitude,  while  in  diffusive  profiles,  the  decay  rate  depends  on  the 
molecular  mass  of  the  constituent  (heavier  ones  decay  more  rapidly).  The  mixing  profiles  dominate 
below  105  km,  and  the  diffusive  profiles  dominate  above,  with  a  smooth  transition  between  the  two 
regions. 

The  MSIS-83  empirical  model  {Hedin,  1983]  was  based  on  data  from  seven  satellites  and  numerous 
rocket  probes,  as  well  as  five  ground  based  incoherent  scatter  stations.  This  model  generated 
predictions  of  temperature,  along  with  densities  for  Nj,  O2,  O,  He,  Ar,  and  H.  With  the  availability 
of  temperature  and  composition  measurements  from  the  Dynamics  Explorer  B  (DE-B),  a  revised 
version  of  the  model  was  created,  MSIS-86  [He(Un,  1987].  The  DE-B  data  provided  good  coverage 
of  the  polar  regions  just  after  the  peak  of  solar  cycle  21,  and  suggested  several  refinements  in  the 
description  of  polar  region  morphology.  Based  on  the  DE-B  data,  plus  data  from  Atmospheric 
Explorers  C,  D,  and  E,  atomic  nitrogen  was  added  to  the  list  of  species.  Both  MSIS-83  and  MSIS-86 
had  a  low  altitude  limit  near  the  mesopause  (85  km). 

The  MSIS-90  thermospheric  model  [Hedin,  1991]  revises  MSIS-86  in  the  lower  thermosphere  and 
extends  the  model  into  the  middle  and  lower  atmosphere,  providing  a  single  analytic  model  of 
temperature  and  density  profiles.  This  extension  through  the  lower  atmosphere  is  based  mainly  on 
tabulations  from  the  Middle  Atmosphere  Program  (MAP)  Handbook  16  by  Barnett  and  Comey 
[1985].  Rocket  and  incoherent  scatter  data  in  the  upper  mesosphere  and  lower  thermosphere  have 
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also  been  incorporated  in  the  revision.  Latitude,  annual,  semiannual,  and  simplified  local  time  and 
longitude  variations  are  modeled  using  low-order  spherical  harmonics  and  Fourier  series. 

The  sample  MSIS-90  output  shown  in  Table  17  was  generated  for  1994,  with  a  geodetic  latitude  of 
40°,  and  a  longitude  of  0°.  Two  solar  flux  and  two  geomagnetic  activity  cases  are  shown  for  each 
altitude/day/local  time  combination.  Note  that  in  addition  to  total  mass  density,  the  model  can  also 
compute  the  following: 


•  Number  densities  (per  cm^)  for  Ar,  H,  He,  N,  Nj,  O,  and  Oj 

•  Exospheric  temperature 

•  Temperature  at  given  altitude 


A  test  on  a  66  MHz  486  PC  making  10,000  calls  to  the  MSIS-90  routine  required  49  seconds;  the 
same  test  ran  in  40  seconds  on  a  40  MHz  Sparc  IPX  workstation.  A  capability  of  using  an 
extensive  ap/Fjoj  database  also  exists;  this  is  useful  when  making  "historical"  model  runs. 


Table  17.  Sample  MSIS-90  Output 


Alt 

Day 

LT 

ap  =  7 
Fio.7  ~  150 

(Kp  =  2) 

Fio.7  —  220 

ap  =  48 
Fio.7  ~  150 

(Kp  =  5) 

Fio.7  =  220 

0 

1 

0 

1.247E-03 

1.237E-03 

1.244E-03 

1.243E-03 

0 

1 

12 

1.237E-03 

1.237E-03 

1.243E-03 

1.242E-03 

0 

90 

0 

1.214E-03 

1.215E-03 

1.224E-03 

1.223E-03 

0 

90 

12 

1.215E-03 

1.215E-03 

1.224E-03 

1.223E-03 

0 

180 

0 

1.199E-03 

1.200E-03 

1.210E-03 

1.209E-03 

0 

180 

12 

1.201E-03 

1.201E-03 

1.210E-03 

1.209E-03 

300 

1 

0 

2.014E-14 

3.789E-14 

2.511E-14 

4.565E-14 

300 

1 

12 

2.833E-14 

5.236E-14 

3.559E-14 

6.364E-14 

300 

90 

0 

2.367E-14 

4.213E-14 

2.871E-14 

4.919E-14 

300 

90 

12 

3.544E-14 

6.180E-14 

4.257E-14 

7.176E-14 

300 

180 

0 

1.813E-14 

3.301E-14 

2.262E-14 

3.915E-14 

300 

180 

12 

3.021E-14 

5.295E-14 

3. 629E-14 

6.098E-14 

600 

1 

0 

9.373E-17 

3.165E-16 

1.303E-16 

4.231E-16 

600 

1 

12 

2.164E-16 

7.696E-16 

3.056E-16 

1.031E-15 

600 

90 

0 

1.279E-16 

4.168E-16 

1.854E-16 

5.576E-16 

600 

90 

12 

3.709E-16 

1.226E-15 

5.180E-16 

1.584E-15 

600 

180 

0 

7.228E-17 

2.560E-16 

1.137E-16 

3.560E-16 

600 

180 

12 

2.809E-16 

9. 687E-16 

3.991E-16 

1.243E-15 
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6.2  Groves/MSIS  Atmospheric  Model 

The  Groves/MSIS  model  is  a  hybrid,  similar  to  the  AFRA-86  model.  The  model  is  defined  by 
a  Global  Reference  Atmosphere  [Groves,  1985]  for  altitudes  from  18  to  70  km,  and  by  the  MSIS- 
83  thermospheric  model  above  130  km.  A  transitional  model  is  used  in  the  region  between  70  and 
130  km. 

The  Global  Reference  Atmosphere  is  derived  from  monthly  means  of  zonal  mean  temperature, 
pressure,  density,  number  density,  and  pressure  scale  height.  These  means  have  been  tabulated 
at  heights  fi-om  18  to  80  km,  and  latitudes  from  80°S  to  80‘’N,  with  a  latitude  interval  of  10°. 
Monthly  mean  longitudinal  variations  of  temperature,  pressure,  and  density  have  also  been 
tabulated  for  specific  months  in  each  hemisphere.  The  zonal  means  are  derived  from  tabulations 
of  temperature  and  geopotential  height  based  on  the  following: 

•  Nimbus  5  SCR  and  Nimbus  6  PMR. 

•  A  southern  hemisphere  reference  atmosphere  prepared  for  the  1984  COSPAR  meeting  in 
Graz,  Austria,  based  on  rocketsonde  data. 

•  Two  earlier  northern  hemisphere  rocket-based  reference  atmospheres,  CIRA-72  and  the  Air 
Force  Reference  Atmosphere  1978  (AFRA-78). 

The  longitudinal  variations  are  derived  solely  from  satellite  based  tabulations. 

The  transitional  model  [Groves,  1987]  was  designed  to  maintain  continuity  between  the  Global 
Reference  Atmosphere  and  MSIS-83  in  the  second  derivatives  with  respect  to  height  of 
temperamre,  pressure,  density,  and  constituent  gas  concentrations.  In  addition,  a  best  fit  is 
obtained  to  available  temperature  data,  taken  at  5  km  steps  from  75  km  to  125  km.  This 
temperature  data  consisted  of  that  used  for  construction  of  CIRA-72,  and  rocket  and  incoherent 
scatter  data  presented  at  the  1984  COSPAR  meeting  in  Graz,  Austria.  The  transitional  model  was 
also  designed  to  reproduce  the  required  ratio  of  N2  pressure  at  70  km  to  that  at  130  km,  on 
integration  of  appropriate  physical  equations. 


6.3  Jacchia  Thermospheric  Models 

The  Jacchia  models  [Jacchia,  1964;  1970;  1971;  1977]  generate  temperatures  and  densities  using 
analytic  vertical  temperature  profiles.  These  temperature  profiles  depend  on  geographic 
coordinates,  local  time,  day  of  year,  and  geophysical  conditions  (solar  flux  and  geomagnetic 
activity).  Exospheric  temperature  is  defined  by  these  input  parameters;  the  vertical  profile  is  then 
determined  from  the  exospheric  temperature,  with  a  geomagnetic  activity  dependent  correction 
in  the  1977  model. 

The  1964  model  generates  profiles  starting  at  120  km,  and  is  based  primarily  on  satellite  drag 
measurements.  Analytic  temperature  profiles  are  generated  based  on  fixed  boundary  conditions 
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at  120  km.  These  boundary  conditions  constitute  a  major  drawback,  since  temperature  and  density 
undergo  considerable  variations  at  120  km.  The  1970  version  is  based  on  boundary  conditions 
at  90  km,  and  contains  a  higher  ratio  of  atomic-oxygen  to  molecular-oxygen  density.  Mixing  of 
the  various  constituents  (nitrogen,  oxygen,  argon,  helium,  and  hydrogen)  is  assumed  up  to  105 
km,  with  diffusion  above  this  height.  Moving  the  boundary  to  90  km  also  required  a  change  in 
the  analytic  expression  used  to  generate  the  temperature  profiles. 

The  1970  model  was  approved  for  inclusion  in  the  COSPAR  International  Reference  Atmosphere 
(CIRA),  subject  to  some  suggested  revisions.  The  molecular  oxygen  concentration  was  decreased 
further  and  that  of  atomic  oxygen  was  increased,  bringing  the  composition  at  150  km  in  line  with 
the  results  of  a  survey  of  experimental  data.  The  transition  from  a  condition  of  mixing  of  the 
various  constituents  to  a  state  of  diffusion  was  shifted  to  100  km.  This  revised  version  of  the 
Jacchia  model  was  incorporated  in  CIRA  [COSPAR,  1972],  and  is  known  separately  as  the  Jacchia 
1971  model. 

The  1977  Jacchia  model  (J77)  is  based  on  satellite  drag  measurements  taken  from  1960  to  1975 
above  200  km.  In  addition,  the  1977  model  incorporates  mass  spectrometer  measurements  from 
satellites  OGO  6  and  ESRO  4.  An  attempt  is  made  to  reproduce  the  results  from  OGO  6  for  the 
relative  concentrations  of  N2  and  O  at  450  km,  while  keeping  the  total  density  profiles  anchored 
to  satellite  drag.  A  similar  approach  was  taken  in  the  construction  of  the  high  altitude  portion  of 
the  1976  U.S.  Standard  Atmosphere  [COESA,  1976].  Consequently,  the  J77  profQes  at  1000  K 
are  very  similar  to  the  U.S.  Standard  profiles,  which  are  defined  for  an  exospheric  temperature 
of  1000  K. 

Note  that  although  the  data  coverage  in  altitude  and  composition  is  not  as  extensive  as  for  the 
MSIS  models,  the  solar  activity  reached  a  higher  level  during  the  data  collection  period  than  that 
reached  during  the  cycle  covered  by  most  of  the  MSIS  measurements.  Therefore,  the  solar 
activity  coverage  of  the  Jacchia  data  is  more  extensive.  A  lack  of  flexibility  in  the  vertical 
temperature  profiles  results  in  poor  representation  of  variations  with  local  time  and  geomagnetic 
activity.  The  Jacchia  models  lack  a  semidiurnal  local  time  variation  (two  oscillations  per  day), 
normally  seen  at  low  altimdes  but  not  at  high  altitudes.  This  is  a  result  of  using  temperature 
profiles  depending  only  on  exospheric  temperature  (that  is,  high  altitude  data). 

The  sample  Jacchia  77  output  shown  in  Table  18  was  generated  for  1994,  with  a  latitude  of  40° 
and  a  longitude  of  0°.  Two  solar  flux  and  two  geomagnetic  activity  cases  are  shown  for  each 
altimde/day/local  time  combination.  A  test  on  a  40  MHz  Sparc  IPX  workstation  making  10,000 
calls  to  the  Jacchia  77  routine  required  21  seconds. 
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Table  18.  Sample  Jacchia  77  Output 


tJ 

II 

(Kp  =  2) 

ap  =  48 

(Kp  =  5) 

Day 

LT  . 

Alt 

Fio.7  =  150 

Fio.7  =  220 

Fio.7  =  150 

Fio.7  =  220 

1 

0 

100 

5.832E-10 

5. 821E-10 

6.340E-10 

6.328E-10 

1 

12 

100 

5.821E-10 

5. 811E-10 

6.328E-10 

6.317E-10 

90 

0 

100 

6.278E-10 

6.266E-10 

6.833E-10 

6.820E-10 

90 

12 

100 

6.272E-10 

6.260E-10 

6.826E-10 

6.814E-10 

180 

0 

100 

5.528E-10 

5.518E-10 

6.021E-10 

6.010E-10 

180 

12 

100 

5.519E-10 

5.509E-10 

6.012E-10 

6.001E-10 

1 

0 

300 

1.724E-14 

2.933E-14 

1.833E-14 

3.093E-14 

1 

12 

300 

2.515E-14 

3.998E-14 

2.662E-14 

4.208E-14 

90 

0 

300 

2.495E-14 

4.096E-14 

2.695E-14 

4.404E-14 

90 

12 

300 

3.401E-14 

5.285E-14 

3.668E-14 

5.676E-14 

180 

0 

300 

1.938E-14 

3. 123E-14 

2.124E-14 

3.406E-14 

180 

12 

300 

2.544E-14 

3.903E-14 

2.784E-14 

4.250E-14 

1 

0 

600 

7.437E-17 

2.405E-16 

8.305E-17 

2.587E-16 

1 

12 

,  600 

1.675E-16 

5.043E-16 

1.850E-16 

5.307E-16 

90 

0 

600 

1 . 568E-16 

5.086E-16 

1.771E-16 

5.389E-16 

90 

12 

600 

3.253E-16 

9.455E-16 

3.564E-16 

9.811E-16 

180 

0 

600 

1.165E-16 

3.548E-16 

1.289E-16 

3.713E-16 

180 

12 

600 

2.204E-16 

6.075E-16 

2.376E-16 

6.259E-16 

6.4  Air  Force  Reference  Atmosphere  - 1986 

The  Air  Force  Reference  Atmosphere  (AFRA)  [Bass,  et  al.,  1987]  specifies  temperature, 
pressure,  and  total  mass  density  between  80  and  200  km.  The  AFRA  model  was  derived  by 
merging  the  MSIS-83  model  (defined  above  85  km)  with  a  model  of  Forbes  [1985],  which  is 
based  on  low  altitude  rocket  and  incoherent  scatter  measurements  of  temperature.  The  Forbes 
model  is  used  in  the  lower  portion  of  the  AFRA  altitude  range,  while  MSIS-83  is  used  in  the 
upper  portion.  A  transition  region  provides  a  continuous  and  smooth  match  between  the  two 
models. 

The  Forbes  model  specifies  monthly  mean  profiles  of  total  mass  density,  temperature,  and 
pressure  between  65  and  120  km  in  tabular  form.  Temperature  profiles  were  generated  first  by 
least  squares  fits  to  the  data.  The  pressure  profiles  were  obtained  as  solutions  to  the  barometric 
equation,  using  the  pressure  at  68  km  given  by  the  AFRA-78  [Cole,  et  al,  1978]  as  a  lower 
boundary  condition.  The  density  is  then  derived  from  the  perfect  gas  law. 

The  Forbes  model  depends  only  on  altitude,  latitude,  and  month,  while  MSIS-83  has  additional 
dependencies  on  local  time,  solar  activity,  and  geomagnetic  activity.  To  determine  the  boundaries 
of  the  transition  region,  the  additional  dependencies  were  removed  from  MSIS-83 .  Comparisons 
of  the  two  models  showed  a  distinct  disagreement  between  the  two  models,  which  worsened  going 
fi:om  low  to  mid  latitude.  Therefore,  it  was  concluded  that  the  two  models  could  not  be  brought 
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into  agreement  by  a  minor  adjustment  without  compromising  their  validity  in  their  separate 
regions  of  application. 

The  model  definition  that  was  adopted  uses  the  Forbes  model  below  104  km  and  MSIS-83  above 
120  km.  A  connecting  model  was  defined  for  the  region  between  104  and  120  km.  This 
connecting  model  computes  the  density,  temperature,  and  molecular  weight  from  cubic 
polynomials  matching  the  Forbes  model  at  102  and  104  km,  and  matching  MSIS-83  at  120  and 
122  km.  The  pressure  is  then  computed  from  the  perfect  gas  law. 


6.5  U.S.  Standard  Atmosphere  1976 

The  U.S.  Standard  Atmosphere,  1976  (US76)  was  adopted  by  the  United  States  Committee  on 
Extension  to  the  Standard  Atmosphere  (COESA)  in  February,  1975  [1976].  This  version  is  the 
same  as  the  1962  Standard  Atmosphere  below  50  km,  but  substitutes  newer  values  at  higher 
altitudes.  COESA  based  the  1976  revision  on  the  World  Meteorological  Organization’s  (WMO) 
definition  of  a  standard  atmosphere. 

US76  is  a  steady-state  representation  of  mean  annual  conditions  of  the  atmosphere  at  a  latitude  of 
45°N,  defined  for  periods  of  moderate  solar  activity.  The  model  is  defined  by  sea-level 
temperature  and  pressure,  and  a  temperature-height  profile  to  1000  km.  The  lower  32  km  of  the 
profile  is  identical  to  the  International  Civil  Aviation  Organization  (ICAO)  Standard  Atmosphere. 
The  portion  between  32  and  55  km  is  based  on  radiosonde  and  meteorological  rocketsonde 
observations.  From  55  to  86  km,  the  model  is  based  primarily  on  measurements  from  grenade, 
pitot  static  mbe,  and  falling  sphere  experiments.  Above  86  km,  four  functions  are  used  that  yield 
a  continuous  first  derivative  with  respect  to  altitude  over  the  entire  region.  These  functions  are 
based  on  various  kinds  of  observational  data  [Jursa,  1985]. 

The  basis  for  the  1976  revision  is  a  much  greater  collection  of  experimental  data,  some  over  parts 
of  the  solar  cycle  not  available  for  the  1962  version.  Statistics  compiled  during  preparation  of  the 
1976  version  showed  that  densities  were  about  10  percent  lower  in  the  70  to  80  km  region,  and 
10  percent  higher  in  the  90  km  region  than  those  in  the  1962  Standard.  The  1976  version  uses 
an  exospheric  temperature  of  1000  K,  compared  to  1500  K  in  the  1962  model. 

Presently,  the  US76  code  will  generate  temperamre,  pressure,  molecular  weight,  and  total  mass 
density  up  to  500  km.  A  test  making  1 ,000,000  calls  to  the  US76  routine  required  24  seconds  on 
a  40  MHz  Sparc  IPX  workstation.  A  typical  US76  density  profile  is  shown  in  Table  19. 
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Table  19.  Sample  US76  Output 


Alt 

Density 

Temp 

Pressure 

Mole  Wt 

(km) 

(g/cm^) 

(°  K) 

(dynes/cm^) 

(g/Mole) 

100 

5.604E-10 

195.08 

3.200E-01 

28.40 

150 

2.076E-12 

634.39 

4.543E-03 

24.10 

200 

2.541E-13 

854.56 

8.476E-04 

21.30 

250 

6.073E-14 

941.33 

2.477E-04 

19.19 

300 

1.916E-14 

976.01 

8.769E-05 

17.73 

350 

7.014E-15 

990.06 

3.449E-05 

16.74 

400 

2.803E-15 

995.83 

1.452E-05 

15.98 

450 

1.184E-15 

998.22 

6.443E-06 

15.25 

The  pressure  values  produced  by  the  model  code  are  calculated  from  the  ideal  gas  law  instead  of 
the  original  tables;  small  differences  due  to  round-off  errors  may  occur. 


6.6  Global  Reference  Atmospheric  Model  1988 

The  Global  Reference  Atmospheric  Model  1988  (GRAM88)  is  a  hybrid  consisting  of  three  distinct 
zones  and  two  transition  regions  [Justus,  et  al,  1980] .  This  model  was  developed  by  Georgia 
Tech,  under  NASA  sponsorship.  GRAM88  is  valid  over  an  altitude  range  of  0  to  700  km,  and 
contains  variations  with  latitude,  longitude,  and  month.  Table  20  summarizes  the  usage  of  the 
Jacchia  1970  and  Groves  1971  models,  and  a  4-D  atmospheric  model  developed  for  NASA  by 
Allied  Research  Associates. 


115 

90 

30 

25 

Surface 


Table  20.  GRAM88  Model  Composition 
Jacchia  1970 

Transition  between  Jacchia  and  Groves  models 
Groves  1971 

Interpolation  between  4-D  data  and  Groves  model 
NASA  4-D  Model 


The  4-D  atmospheric  model  is  based  on  empirically  determined  atmospheric  parameter  profiles 
taken  at  a  large  array  of  locations.  This  data  consists  of  mean  monthly  and  daily  variance  profiles 
of  pressure,  density,  and  temperature  at  1  km  intervals  from  the  surface  to  25  km,  for  the  entire 
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globe.  These  data  are  interpolated  over  latitude  and  longitude  to  find  the  necessary  values.  In  the 
transitional  region  from  25  to  30  km,  an  interpolation  over  height  is  performed  between  the  NASA 
4-D  model  and  the  Groves  model. 

The  Groves  1971  empirical  model  combines  many  observations  from  a  wide  variety  of  longitudes, 
using  observational  results  over  a  period  of  approximately  6  years.  Latitude  coverage  of  the 
Groves  model  is  from  the  equator  to  70®  -  80°.  Southern  hemisphere  data  is  taken  as  northern 
hemisphere  data  with  a  six  month  change  of  date.  The  Groves  model  data  contains  no  longimde 
variation;  rather,  the  data  is  modified  by  longitude,  latitude,  and  height  dependent  stationary 
perturbations. 

A  smooth  transition  between  the  Groves  and  Jacchia  models  is  obtained  using  a  "fairing" 
technique,  described  by  Justus,  et  al,  [1974].  This  fairing  is  performed  only  at  95,  100,  105, 
and  110  km,  which  are  the  altitudes  for  which  there  are  Groves  values.  The  remaining  locations 
in  this  transition  region  are  then  filled  in  using  linear  interpolation. 

The  sample  GRAM-88  ouq)ut  shown  in  Table  21  was  generated  for  1994,  with  a  geodetic  latitude 
of  40°,  and  a  longimde  of  0°.  Two  solar  flux  and  two  geomagnetic  activity  cases  are  shown  for 
each  altimde/day/local  time  combination.  A  test  run  calculating  24  profiles,  each  containing  three 
altimde  points,  required  288  seconds  on  a  40  MHz  Sparc  IPX  workstation. 


Table  21.  Sample  GRAM-88  Output 


ap  =  7 

Alt  Day  LT  Fio.7  ~  150 


0 

1 

0 

0 

1 

12 

0 

90 

0 

0 

90 

12 

0 

180 

0 

0 

180 

12 

300 

1 

0 

300 

1 

12 

300 

90 

0 

300 

90 

12 

300 

180 

0 

300 

180 

12 

600 

1 

0 

600 

1 

12 

600 

90 

0 

600 

90 

12 

600 

180 

0 

600 

180 

12 

1.246E-03 
1.246E-03 
1.235E-03 
1.235E-03 
1.185E-03 
1.185E-03 
1. 815E-14 
2.747E-14 
2.312E-14 
3.469E-14 
1. 902E-14 
3. 008E-14 
7.353E-17 
1. 964E-16 
1.296E-16 
3.513E-16 
8.192E-17 
2.459E-16 


(Kp  -  2) 
Fio.7  =  220 

1.246E-03 
1.246E-03 
1.235E-03 
1.235E-03 
1.185E-03 
1.185E-03 
2.396E-14 
3.362E-14 
2.918E-14 
4.081E-14 
2 . 488E-14 
3. 624E-14 
1.410E-16 
3.243E-16 
2.280E-16 
5 .317E-16 
1.546E-16 
3 . 924E-16 


ap  =  48 
Fio.7  =  150 

1.246E-03 
1.246E-03 
1.235E-03 
1.235E-03 
1.185E-03 
1.185E-03 
3.216E-14 
4.431E-14 
3, 960E-14 
5.351E-14 
3.286E-14 
4 .700E-14 
2.901E-16 
6.588E-16 
4.922E-16 
1.092E-15 
3 . 065E-16 
7 .700E-16 


(Kp  =  5) 
Fio.7  =  220 

1.246E-03 
1.246E-03 
1.235E-03 
1.235E-03 
1. 185E-03 
1.185E-03 
3.831E-14 
5.011E-14 
4.559E-14 
5.881E-14 
3.900E-14 
5.267E-14 
4.517E-16 
9.135E-16 
7.101E-16 
1. 420E-15 
4.733E-16 
1.046E-15 
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7.  ORBITAL  PREDICTION  MODELS 


Orbital  prediction  models  provide  the  basis  for  the  analysis  of  space-based  observations  of 
astronomical  interest  as  well  as  Earth  and  near  space  environments.  Where  high  precision  is  not 
required,  orbital  prediction  based  upon  NORAD  orbital  elements,  provided  by  SPACECOM  and 
CSTC  (SCF),  are  suitable  for  most  purposes.  The  first  of  the  orbital  prediction  models  described 
here  (LOKLANGL)  has  been  developed  for  use  at  Phillips  Laboratory,  and  has  been  customized 
for  use  with  a  number  of  PL  projects,  in  particular  those  involving  Earth’s  magnetic  field 
(ionosphere  and  magnetosphere).  LOKANGL  uses  an  analytic  perturbation  model.  The  second 
model  is  the  standard  Spacetrack  package  (SGP,  SGP4,  etc.)  which  implements  five  related 
analytic  perturbation  models.  A  comparison  of  LOKANGL  with  the  Spacetrack  models  has  been 
performed  for  several  satellites  pertinent  to  PL  applications,  with  the  finding  that  LOKANGL 
satisfactorily  matches  the  recommended  Spacetrack  model  (as  determined  by  the  satellite’s  orbital 
parameters). 

To  obtain  higher  precision  than  the  LOKANGL  and  the  Spacetrack  models,  it  is  usually  necessary 
to  implement  a  numerical  integration  scheme  using  a  gravitational  model  (a  spherical  harmonic 
representation  of  Earth’s  gravitational  field)  and,  for  near  earth  satellites,  an  atmospheric  density 
model,  to  compute  the  effect  of  drag.  When  tracking  data  are  available,  orbital  parameters  may 
be  optimized  using  statistical  weighting  techniques. 

The  third  model,  ASAP,  developed  by  JPL,  provides  this  capability.  At  PL  an  adaptation  of  the 
ASAP  program  has  been  used  to  investigate  lifetimes  of  satellites  whose  orbital  decay  is  governed 
mostly  by  atmospheric  drag.  The  latter  is  highly  responsive  to  spatial  and  dynamical  variations 
in  the  neutral  density.  Some  of  the  density  models  used  in  these  studies  attempt  to  model 
variations  due  to  solar-terrestrial  effects  due  to  extreme  ultra  violet  (EUV)  radiation  and 
geomagnetic  field  heating  due  to  charged  particle  precipitation. 

Generally,  the  orbital  elements  provided  by  SPACECOM  and  CSTC  are  limited  to  the  accuracy 
of  the  Spacetrack  models.  To  obtain  higher  precision  than  these  three  models  permit  requires  an 
orbital  analysis  package,  which  includes  not  only  the  capability  to  predict  the  position  of  a 
satellite,  but  also  to  compare  predictions  with  actual  observations,  and  perform  differential 
correction,  with  either  availability  of  such  observations,  or  the  capability  and  resources  to  make 
such  observations. 


7.1  LOKANGL  Orbital  Prediction  Program 

LOKANGL  [Hein,  et  al. ,  1991]  is  an  orbit  analysis  and  prediction  program  that  has  served  broad 
requirements  at  Phillips  Laboratory  for  many  years.  This  program  predicts  the  position  and  other 
orbital  parameters  of  an  Earth  satellite  using  either  NORAD  or  CSTC  (SCF)  orbital  elements. 
LOKANGL  also  provides  estimates  of  several  important  geophysical  parameters  at  the  location 
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of  the  satellite,  such  as  L-shell  value  and  dipole  or  corrected  magnetic  latitude  and  local  time.  In 
addition  to  the  standard  orbit  prediction  mode,  LOKANGL  employs  an  interpolation  mode  when 
generating  an  orbit  that  is  bounded  in  time  by  a  pair  of  successive  element  sets.  Augmentation 
and  adaptation  of  LOKANGL  continues  in  response  to  various  needs. 

In  the  ideal  case  of  two  body  motion,  the  orbit  of  two  spherically  symmetric  masses  bound  by 
mutual  gravitational  attraction  is  an  ellipse.  However,  for  Earth  satellites,  deviations  of  Earth’s 
mass  distribution  from  spherical  symmetry  must  be  considered.  LOKANGL  models  this  deviation 
using  an  analytic  perturbation  technique,  which  takes  into  account  the  effect  of  the  equatorial 
bulge  and  the  north-south  asymmetry. 

In  general,  precision  orbit  prediction  programs  numerically  integrate  the  equations  of  motion.  By 
contrast,  LOKANGL  uses  an  analytic  perturbation  computation  that  is  much  more  efficient  in 
terms  of  computer  time.  This  technique  takes  into  account  the  second  and  third  zonal  harmonics 
of  the  Earth’s  gravitational  field.  Omission  of  the  high  order  gravitational  perturbation  terms 
results  in  periodic  errors  of  a  few  kilometers  in-track  within  a  single  orbit.  Other  sources  of  error 
include  the  effect  of  atmospheric  drag  and  lunar-solar  perturbations.  A  study  has  been  made  to 
determine  the  effect  of  spacing  of  orbital  elements  on  in-track  errors  [Hein  and  Robinson,  1991], 
and  may  be  consulted  for  further  information. 

The  major  capabilities  of  LOKANGL  are  the  following: 

•  Computation  of  a  satellite’s  position  and  velocity  vectors  at  regularly 
spaced  time  intervals. 

•  Calculation  of  a  satellite’s  mean  elements  at  regularly  spaced  time  intervals. 

•  Determination  of  a  satellite’s  position  relative  to  as  many  as  20  stations 
located  on  the  Earth’s  surface. 

•  Trajectory  interpolation  using  pairs  of  orbital  elements. 

•  Limited  handling  of  orbital  maneuvers  is  available. 

•  Modeling  of  atmospheric  density/drag  and  geopotential  for  evaluating  time 
derivatives  of  elements  based  on  the  method  of  King-Hele. 

•  Although  not  presently  embedded  in  the  program  due  to  infrequent  usage, 
a  capability  of  providing  look  angles  to  satellites  from  moving  stations 
(such  as  airplanes)  is  available. 

•  Calculation  of  solar  illumination  and  shadowing  of  satellites. 

•  Options  for  printing  a  variety  of  output  listings. 

•  A  menu  interface  is  available  with  the  PC  version  of  LOKANGL. 

Table  22  gives  position  and  velocity  vectors  for  Mir  over  a  2-hour  interval  at  10-minute 
increments,  while  Table  23  presents  sub-satellite  data  for  the  same  time  period.  In  a  test  on  a  66 
MHz  486  PC,  calculation  of  sub-satellite  data  at  1-minute  steps  for  a  5-day  period  required  50 
seconds.  Generation  of  position  and  velocity  vectors  over  the  same  period  required  42  seconds. 
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Table  22.  Position  and  Velocity  Vectors 


DATE 

U  TIME 

SATELLITE  16609 

MO  DY 

YR 

HR  m  sc 

X(KM)  Y(KM) 

Z(KM) 

XD0T(KM/SEC) 

YDOT(KM/SEC) 

ZDOTC KM/SEC) 

6 

8 

94 

0  0 

0 

•’3446.0984  2422.8861 

-5314.7982 

-4.575680 

-6.141422 

.169238 

6 

8 

94 

0  10 

0 

-5225.2233  -1521.4308 

-4043.9981 

-1.124580 

-6.499486 

3.904578 

6 

8 

94 

0  20 

0 

-4691.7660  -4791.4141 

-979.2058 

2.836488 

-3.977076 

5.916245 

6 

8 

94 

0  30 

0 

-2075.4085  -5933.4384 

2520.8666 

5.546864 

.318820 

5.296606 

6 

8 

94 

0  40 

0 

1463.5880  -4437.8990 

4897.3146 

5.792170 

4.472043 

2.317365 

6 

8 

94 

0  50 

0 

4352.8755  -971.1934 

5092.8746 

3.465825 

6.635685 

-1.690460 

6 

8 

94 

1  0 

0 

5309.8816  2926.5617 

3021 .4489 

-.399898 

5.853232 

-4.946308 

6 

8 

94 

1  10 

0 

3908.8854  5524.6458 

-395-4133 

-4.088609 

2.471850 

-6.000387 

6 

8 

94 

1  20 

0 

773.1740  5671.7533 

-3636.9303 

-5.958497 

-1.998035 

-4.386461 

6 

8 

94 

1  30 

0 

-2705.4529  3309.2061 

-5265.5130 

-5.190155 

-5.572081 

- .833659 

6 

8 

94 

1  40 

0 

-4989.5734  -514.6431 

-4562.6656 

-2.129308 

-6.683162 

3.087027 

6 

8 

94 

1  50 

0 

-5067.0546  -4109.9947 

-1837.0915 

1.883670 

-4.836818 

5.647690 

6 

8 

94 

2  0 

0 

-2896.5254  -5881.0159 

1705.0441 

5.071887 

-.836330 

5.700232 

Table  23. 

Sub-Satellite  Output 

DATE 

U  TIME 

GEOD 

SATELLITE  16609 

MO  DY  YR  HR 

MN  SC 

REV 

LAT  E  LONG  ALT(KM) 

SZ  ANG  SUNSH 

LUNSH 

LOCT  GMLAT 

GMLONE 

GMLT  L 

-SHELL  B/BO 

6 

8 

94 

0 

0 

0  47460 

-51.775  -111.260  416.810 

95.66  14.4 

6.5 

16.58  -43.23 

326.75 

17.32 

1.944  9.16 

6 

8 

94 

0 

10 

0  47460 

-36.789  -62.423  409.704 

127.22  -17.2 

-29.4 

20.01  -26.26 

7.69 

20.21 

1.344  1.75 

6 

8 

94 

0 

20 

0  47460 

-8.360  -35.561  399.389 

147.51  -37.8 

-53.0 

21.96  .33 

35.05 

22.19 

1.143  1.06 

6 

8 

94 

0 

30 

0  47461 

21.978  -12.949  397.386 

135.06  -25.2 

-27.8 

23.64  27.27 

62.34 

.17 

1.143  1.45 

6 

8 

94 

0 

40 

0  47461 

46.523  22.075  402.155 

104.91  5.1 

8.3 

2.14  44.93 

103.95 

3.10 

1.931  9.44 

6 

8 

94 

0 

50 

0  47461 

48.970  78.738  403.609 

72.10  37.8 

45.3 

6.08  39.55 

154.63 

6.64 

1.977  11.82 

6 

8 

94 

1 

0 

0  47461 

26.634  117.671  400.230 

42.71  67.1 

80.2 

8.84  16.11 

188.10 

9.03 

1.112  1.71 

6 

8 

94 

1 

10 

0  47461 

-3.365  141.022  401.100 

33.37  76.5 

85.8 

10-57-12.34 

212.86 

10.84 

1.055  1.30 

6 

8 

94 

1 

20 

0  47461 

-32.594  166.033  409.882 

55.57  54.3 

52.7 

12.40  -37.87 

243.63 

13.05 

1.783  8.38 

6 

8 

94 

1 

30 

0  47461 

-51.108  -149.444  416.781 

87.36  22.7 

16.1 

15.54  -47.87 

293.52 

16.53 

2.506  23.34 

6 

8 

94 

1 

40 

0  47461 

-42.470  -95.329  412.296 

119.56  -9.5 

-20.0 

19.31  -32.65 

338.88 

19.71 

1.470  3.15 

6 

8 

94 

1 

50 

0  47461 

-15.821  -64.679  401.466 

144.53  -34.7 

-50.0 

21.52  -5.25 

6.15 

21.69 

1.143  1.04 

6 

8 

94 

2 

0 

0  47462 

14.667  -42.453  396.924 

140.75  -30.9 

-36.2 

23.17  23.93 

30.43 

23.46 

1.211  1.66 

7.2  Spacetrack  Models  for  Propagation  of  NORAD  Element  Sets 

The  SPACETRACK  code  implements  five  different  analytic  perturbation  schemes  compatible  with 
the  NORAD  satellite  elements  distributed  by  SPACECOM.  The  NORAD  elements  are  mean 
elements  which  differ  significantly  from  osculating  (instantaneous  position  and  velocity,  or 
equivalent  Keplerean  elements)  by  modeling  periodic  perturbations  due  to  deviations  of  Earth’s 
gravitational  field  from  spherical  symmetry,  and  the  effect  of  Earth’s  atmosphere  upon  satellite 
orbits.  There  is  no  "standard"  definition  of  mean  elements;  each  particular  definition  is 
characterized  by  the  perturbations  which  are  taken  into  account  (modeled)  and  those  which  are 
ignored.  The  SPACETRACK  codes  enable  users  to  make  Position  and  Velocity  predictions  using 
the  NORAD  elements  by  using  perturbation  models  which  are  compatible  with  the  NORAD  mean 
elements.  The  SPACETRACK  code  is  described  in  the  NORAD  document  "SPACETRACK 
REPORT  #3  -  Models  for  Propagation  of  NORAD  Element  Sets",  December  1980. 


According  to  the  SPACETRACK  report,  the  models  recommended  for  element  prediction  are  the 
SGP4  for  near-earth  satellites  (orbital  periods  <  225  minutes)  and  SDP4  for  longer  period 
satellites.  The  latter  model  takes  into  account  certain  "resonance"  terms  in  Earth’s  gravitational 
field  for  synchronous  satellites,  and  for  12-hour  period  satellites,  as  well  as  perturbations  due  to 
the  sun  and  moon.  The  SGP  model  is  an  older  near-earth  perturbation  model  that  has  been 
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replaced  by  SGP4.  In  addition  to  the  SGP,  SGP4  and  SDP4  models,  two  additional  models  are 
available  in  SPACETRACK,  SGP8  and  SDP8,  improved  versions  of  SGP4  and  SDP4,  which, 
according  to  the  document,  have  been  proposed  as  replacements,  but  have  not  been  officially 
adopted. 

Program  input  consists  of  start  and  stop  time,  time  interval  in  minutes  (either  forwards  or 
backwards),  followed  by  a  two  element  set.  Output  is  time  (relative  to  epoch),  and  position 
vectors  and  velocity  vectors  (ECI)  in  [km],  [km/sec]. 

Dr.  Don  Larson,  SAIC  has  provided  the  follovang  details; 

NORAD  elements  are  produced  with  SGP4,  and  are  converted  to  the  composite  SGP/SGP4  2-line 
format  to  be  fully  compatible  with  the  many  SGP  users  (using  the  algorithms  which  have  been 
supplied  to  thousands  of  requestors  down  through  the  years).  The  epoch  time,  eccentricity, 
incWtion,  right-ascension,  argument  of  perigee,  and  mean  anomaly  terms  are  identical  between  SGP 
and  SGP4.  The  SGP4  drag  term  (B-star)  is  included,  as  are  the  SGP  drag  terms  (n-dot/2  and 
n-dot-dot/6).  The  mean  motion  listed  is  the  SGP  mean-mean-motion  (that  is,  Kozai).  It  may  be 
easily  converted  back  to  the  SGP4  (that  is,  Brouwer)  mean  motion  by  removing  the  J2  (Kozai)  term, 
via  the  algorithms  alluded  to  above. 

The  earth  model  and  physical  constants  are  WGS-72,  and  the  coordinate  frame  is:  true  equator  and 
mean  equinox  of  epoch  (as  listed  on  each  element  set),  using  the  FK5  mean  of  J2000  time  and 
reference  frame. 


7.3  ASAP  -  ArtiHcial  Satellite  Analysis  Program 

ASAP  is  a  satellite  prediction  program  using  precision  numerical  integration,  developed  by  JPL 
for  orbiting  planetary  spacecraft.  It  uses  an  eighth  order  Runge-Kutta-Fehlberg  numerical 
integration  routine  with  automatic  step  size  control  for  the  integration  of  the  equations  of  motion 
to  implement  Cowell’s  method  of  special  perturbations. 

ASAP  allows  for  the  use  of  gravitational  potential  spherical  harmonic  expansions  up  to  40th  order 
and  degree,  and  optionally  takes  into  account  lunar  and  solar  gravity,  drag,  solar  radiation 
pressure  (including  planet  shadow  entry  and  exit).  In  addition  to  providing  the  US  1976 
atmospheric  density  model  for  the  Earth,  it  also  allows  for  an  exponential  air  density  model  for 
any  planet.  ASAP  has  been  modified  at  PL  to  utilize  dynamical  atmospheric  density  models  such 
as  Jacchia  70  and  MSISE  90. 

In  addition  to  the  standard  inputs  (epoch  [Ephemeris  Time],  and  orbital  elements  of  the  spacecraft 
in  either  osculating  or  mean  elements),  option  flags,  and  also  model  parameters,  such  as 
coefficients  for  spherical  harmonics  for  gravitational  potential  fields,  or  geophysical  indices  for 
atmospheric  density,  must  be  provided  by  the  user  when  options  are  selected. 
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